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Introduction 

This project was undertaken to examine the approach to steady state in collisional mountain 
belts. Although the primary thrust of this grant was to look at larger collisional mountain belts, 
such as the Himalaya, the Tien Shan, and Southern Alps, we began by looking at smaller 
structures represented by growing and propagating folds. Like ranges that are evolving toward a 
topographic steady state, these folds undergo a series of morphologic changes as they are 
progressively uplifted and eroded. We wanted to document the nature of these changes and to try 
to discern some of the underlying controls on them. We initially focused on the Wheeler Ridge 
anticline in southern California. 

Subsequently, we progressed to looking at the topographic development and the effects of 
differential uplift and glaciation on the Kyrgyz Range in the northern Tien Shan. This range is 
unusual inasmuch as it is transformed along its length from a simple uplift with a largely 
preserved Mesozoic erosion surface arching across it to a highly dissected and heavily glaciated 
uplift in the region where uplift has been sustained at higher rates over longer intervals. 

In efforts to understand the distribution of erosion rates at ltf-lO 5 year time scales, 
cosmogenic radionuclide (CRN) concentrations have been gaining increasingly widespread 
usage (Brown et al., 1995; Riebe et al., 2004; Riebe et al., 2001 ; Vance et al., 2003). Most 
studies to date, however, have been conducted in slowly eroding ranges. In rapidly eroding 
mountains where landslides deliver most of the sediments to the rivers, we hypothesized that 
CRN concentrations could be highly perturbed by the stochastic processes of landsliding. 
Therefore, we undertook the development of a numerical model that simulated the effects of both 
landsliding and grain-by-grain attrition within fluvial catchments. This modeling effort has 
shown the effects of catchment size and erosion rate on CRN concentrations and allows a 
prediction of where to sample to obtain the optimal erosion rate estimates using CRN techniques. 

Finally, we developed computational techniques to operate on DEMs to extract useful 
information that would enable quantification of climate-erosion interactions. In particular, we 
worked on rapid techniques to define catchments of any given range of sizes, to extract channel 
gradients, to combine precipitation information to calculate discharge, and to utilize various 
stream-power models to determine the erosional energy within any given catchment within a 
transect. 

Below, we briefly describe results from Wheeler Ridge, the Kyrgyz Range, the Nepal 
Himalaya, and our numerical modeling. 

Background geology of the Wheeler Ridge, California 

The Wheeler Ridge anticline is the outermost structure in a fold-and-thrust belt (Davis 
and Namson, 1986) that uplifts the San Emigdio Mountains at the interface between the northern 
edge of the Transverse Ranges and the southern end of the San Joaquin valley in Southern 
California (Medwedeff, 1992). Wheeler Ridge is the surface expression of the anticline (Hoots, 
1930), and is an east-west trending, east-plunging anticline whose eastern terminus is buried 
beneath the modem depositional surface (Figure 1). Wheeler Ridge is approximately 10 km long, 
up to 3 km wide, and has a maximum relief of about 500 meters (Medwedeff, 1992). Extensive 
subsurface and surface data supports the interpretation that Wheeler Ridge is underlain by a 
north-directed wedge thrust that merges southward into the basal detachment (Medwedeff, 

1992). Because a 7 ka surface has been uplifted above the modem depositional plain, we infer 
that the Wheeler Ridge anticline is an actively growing structure, even though no recent 
earthquakes have occurred on the Wheeler Ridge Thrust. 


The topography of the Wheeler Ridge anticline comprises abandoned and uplifted 
Holocene and late Pleistocene alluvial surfaces. Particularly at the eastern end of the anticline, 
these surfaces are very well preserved; to the west, the surfaces are increasingly dissected 
leaving remnants on interfluves. The forelimb of the anticline lies on the northern side of 
Wheeler Ridge and is generally steeper and more dissected than the southern side. Previous work 
has dated several of the surfaces using radiocarbon dating and comparative soil stratigraphy 
(Keller et al., 1998, 1999). These results indicate that the anticline is propagating eastward at 
about 29 mm/yr, that the rate of fold uplift of the easternmost portion of the fold is 3 mm/yr and 
that folding inititated about 400 ka ago (Keller et al., 1998, 1999; Medwedeff, 1992). 

The availability of subsurface structural data, dated geomorphic surfaces, aerial 
photographs shot periodically since 1930, and a 10-meter TOPSAR digital elevation model 
makes Wheeler Ridge an ideal place to quantify the tectonic geomorphology of an actively 
growing structure and to calibrate the resolution of the DEM. Below we describe the methods 
used to constrain the temporal and spatial development of topography in response to surficial and 
tectonic processes. 


Geomorphology of the Wheeler Ridge Fold 

Because Wheeler Ridge is propagating to the east, there is a progressively larger amount 
of uplift and erosion from east to west. This shows up very clearly in the geomorphology, 
whereby smooth surfaces characterized by few stream valleys characterize the newly uplifted 
surfaces, whereas deep dissection and steep slopes characterize the older parts of the fold. Using 
the 10-m DEM, we can quantify spatial variations in the magnitude of erosion (Figure 2). These 
clearly indicate that, ignoring the water and wind gaps, the western part of the fold has 
experienced far greater erosion. 

In an effort to quantify the character of and controls on that erosion, we have divided the 
fold surface into forelimb and backlimb surfaces and into age-dependent zones (Figure 3) 
according to the ages and criteria set forth by Keller et al. (1998, 1999). This delineates a suite of 
surfaces ranging from <1 Oky to >200 ka along a 5 km length of the fold. There are systematic 
chnages in the slope angles and hypsometry along the fold (Fig. 4). Hillslope angles steepen, but 
they do not vary linearly with uplift. In fact, the diminished difference in slope angles between 
areas Q4 and Q5 suggest that the threshold angle for landsliding is being approached, such that 
the slopes will not steepen much above this level. The hypsometry shows that, in the older parts 
of the fold, the topography becomes more evenly distributed. Early on, the fold has an uplifted 
crest with short, steep limbs, such that most of the topography is concentrated near the crest. 

Over time, the limbs lengthen and both streams and landslides cut through the uplifted surface, 
such that the topography becomes more evenly distributed. 

We can use the present stream gradients in conjunction with the reconstructed, pristine 
surface (Fig. 2), to estimate the amount of erosion in each stream valley. When we do this 
calculation both for the streams alone and for the entire landscape, we find that there is a 
dramatic increase in erosion rate associated with the older part of the fold (Figure 5). The reasons 
for these differences are not obvious. With the exception of one outlier, the rates show a good 
correlation with mean channel curvature (m/n) and with the mean hillslope gradient. When we 
use a stream power law to predict erosion rates, there is also a generally good correlation (Figure 
6). Surprisingly, we note that within any area with consistent uplift, the north-facing channels are 
always less steep than the south-facing ones. The reason for this is unknown at present. We 
suggest that the drier conditions, more sparse vegetation, and thicker soil carbonate leads to 


higher runoff on the southerly slopes, thereby creating higher stream power for a given slope and 
catchment and permitting the same rate of erosion to be sustained on gentler slopes. We are 
pursuing this further to try to understand what is controlling the erosion patterns along this 
growing fold. 

Pre-Steady-State Topography of the Kyrgy z Range, Tien Shan 

Studies of the geomorphology of the Kyrgyz Range were undertaken to elucidate the role 
of glacial and fluvial weathering on the establishment of steady-state topographic relief. The 
Kyrgyz Range is a 250 km-long, east-west trending range within the Tien Shan of Central Asia. 
Drainage basins within the Kyrgyz Range typically traverse 3 to 4 km of topographic relief, and 
the upper reaches of most of these drainages are presently glaciated. The larger drainage systems 
of the Kyrgyz Range drain to the north, across the active, main frontal thrust system which builds 
the range (Figure 1). South-facing catchments are much smaller with less well-developed 
tributary networks. 

An unconformity formed on resistant Paleozoic metamorphic and crystalline rocks forms 
a prominent pre-uplift marker horizon which is to be used in this study to estimate erosional 
denudation of the Kyrgyz Range. Extensive exposures of this unconformity surface on the east 
end of the range indicate that stream systems in this region are less mature than further west. 

The present study of the Kyrgyz Range exploited this along-strike variation in drainage 
development to explore the evolution of a range toward steady-state topography. Also, extensive 
glaciation in the upper elevation reaches of the Kyrgyz Range provided a strong signal of the role 
of glacial processes in the development of drainage networks. 

Overall, we examined the topographic development and the effects of differential uplift 
and glaciation on the Kyrgyz Range in the northern Tien Shan. This range is unusual inasmuch 
as it is transformed along its length from a simple folded uplift with a largely preserved erosion 
surface arching across it to a highly dissected and heavily glaciated uplift in the region where 
uplift has been sustained at higher rates over longer intervals. 

We have defined the differential uplift using extensive fission-track transects that enable 
us to define the changing height of the partial annealing zone as well as the time of the initiation 
of deformation along strike. The fission-track and topographic results are described in a 
manuscript that has been submitted to Tectonics [Edward R. Sobel, Michael Oskin, Douglas 
Burbank, and Alexander Mikolaichuk, Exhumation of basement-cored uplifts: Example of the 
Kyrgyz Range quantified with apatite fission-track thermochronology: Tectonics, in review]. 
These data yield a clear image of outward propagation from the core of the range. In addition, we 
have mapped the erosion surface and characterized the topographic changes (using the SRTM 
DEM) to define the gradual dissection and removal of the erosion surface in the older (Miocene), 
more uplifted parts of the range, as well as changes in slopes, catchments, and channels with 
increasing dissection. Methods to analyze these topographic data have been developed with the 
intention of rapid application to DEMs anywhere. Thus, considerable effort has been made to 
build a set of portable routines to automatically analyze stream networks. 

For example, a striking attribute of the along-strike topographic variation of the Kyrgyz 
Range is the size and geometry of watersheds from east to west. This trend appears to strongly 
reflect the maturity of the drainage system. Thus, it appears that understanding the topographic 
evolution of the Kyrgyz Range requires a 3-dimensional approach to measurement of stream 
systems. To address this issue, a set of flexible, modular routines were developed that extract 
and analyze all tributary streams of a watershed. Results are maintained within a georeferenced 
framework using Arc-Info GIS so that spatial relationships may be compared within a drainage 
system or with other geologic features. To document the extent of glaciation, data for over 500 


glaciers from the Russian Catalog and georeferenced panchromatic 15m ASTER imagery have 
been integrated into this database. 

Our analyses indicate that glaciation plays a strong role in forming mature drainage 
networks in the Kyrgyz Range. Large, mature drainages within the western half of the study area 
contain extensive regions of glaciated uplands. Within the largest drainage basins, these upland 
areas are expanded along the range crest. Interestingly, the highest peaks within the Kyrgyz 
Range occur not at the range crest, but on the ridges between these mature watersheds. It 
appears that vigorous glaciation of the range crest has promoted an increase in the catchment 
area of a few watersheds, leaving wider ridges between the major watersheds that support higher 
peaks. These ridges also display a pronounced mismatch between downcutting of the trunk 
stream and downcutting of the tributaries. Typically, high-elevation, low relief uplands drain into 
steep tributary streams that then join the gently-sloping central trunk stream. Moraines at the 
outlet of these trunk streams indicate complete glaciation during the Pleistocene and it is likely 
that glacial downcutting is the primary cause of the mismatch between trunk streams and the 
inter-watershed ridge areas. 

In sharp contrast to the large, extensively glaciated north-facing drainages of the Kyrgyz 
Range, the south side of the range is characterized by far less mature, smaller, and steeper 
drainages inset into a partially dissected unconformity surface. Many of these streams have 
poorly-developed concave profiles or slope parallel to the unconformity surface. Even where the 
uppermost reaches of these drainages is occupied by a glacier or cirque, the drainage slope 
appears to be mostly unaffected by glaciation. Two possible causes for the contrast between the 
north and south side of the Kyrgyz Range are proposed: 1 . The southeast side of the range has 
been exhumed much more recently than the more well-developed central part of the range. This 
may be caused by propagation of the frontal thrust system to the east. This however does not 
explain the asymmetry of the drainage network in the central part of the range. 2. A pronounced 
rain-shadow affect has restricted precipitation on the south side of the range. The preliminary 
conclusion from these observations is that formation and intergration of a tributary network is a 
critical step in development of a mature stream profile. 

Our efforts to quantify the effects of glaciation on the development of stream networks 
have yielded new insights on how glaciers modify large-scale topography. Some of these results 
have been submitted to Geology for review [Oskin, M.E., and Burbank, D.W., Alpine landscape 
evolution dominated by cirque retreat: see attached manuscript]. In combination with analysis of 
aerial photos and satellite imagery, we have used the SRTM data extensively to quantify the 
extent of glaciations, changes in snowlines, differences between fluvial and glacial catchments, 
and rates of headward retreat by glaciers. 

A break in relief and stream channel concavity probably related to glacial erosion occurs 
at -2600m within several drainages on the north side of the range. This elevation is 1200-1500m 
lower than the present snowline in the Kyrgyz Range, reported at 3800 to 41 00m for numerous 
glaciers within the range, and is consistent with Pleistocene ELA estimates. The south side of 
the range displays a much less well-developed signature of glacial erosion, as described above. 
Subtle shallowing of stream profiles above 3500m elevation indicates enhanced glacial erosion 
in the uppermost reaches of the south-facing streams. This suggests a considerable difference on 
Pleistocene ELA between the north and south side of the ranges, and lends support to the rain- 
shadow hypothesis proposed above. Analysis of the SRTM data shows that initially streams 
form as shallow linear troughs on the erosion surface. These have little concavity and are 
analogous to halfpipes, with rather uniform valley width along their lengths. As glaciation occurs 


in their headwaters, the upper valley profile flattens and the valley head widens. Where the 
erosion surface is well preserved, we can show that the glaciers initially erode headwards 2-4 
times faster than they erode downwards. This surprising result emphasizes the effects of efficient 
erosion processes in the headvvall area and of relative ineffectual vertical erosion by the juvenile 
glaciers. 

Cosmogenic radionuclide erosion rates in rapidly denuding landscapes 

Samples of fluvial sediment for cosmogenic radionuclides (CRNs) have great potential 
because fluvial samples represent a natural, stochastic integration of all upstream concentrations. 
Landslides, however, have the potential to perturb CRN concentrations away from the long-term 
mean because they can incorporate rocks that have received little cosmic ray exposure. In such 
cases, they can produce low CRN concentrations that suggest rapid denudation. Alternatively, 
because landslides are infrequent events, it is likely that even in small valleys with landslides, 
one has not occurred recently. Hence, when collecting fluvial samples for CRN analysis, one 
should ask, “How large a catchment do I need to sample to iron out all of the CRN irregularities 
due to landsliding?” and “How well does a well-integrated downstream sample represent the 
actual mean of total sediment flux?” To address these questions, we built a numerical model that 
utilizes a DEM, calculates CRN concentrations and production throughout the study area, 
scatters landslides across the area, recalculates CRN distributions in both the remaining bedrock 
and in the landslide debris. We examined different rates of grain-by-grain bedrock weathering in 
combination with various rates of landsliding (0.7 to 9.7 mm/yr that span most of the naturally 
occurring rates. The results of this study are reported in the attached paper that is in press at 
EPSL: Nathan Niemi, Mike Oskin, Douglas Burbanky and Arjun Heimsath, Effects of Bedrock 
Landsliding on Cosmogenically Determined Erosion Rates: EPSL in press: see attached. 

Climate-erosion coupling in the central Himalaya 

We utilized both SRTM and DTED data in the analysis of channels and catchments in the 
central Himalaya. Our related research had discovered that rapid erosion rates persist into the 
rain shadow of the Himalaya. So, despite a 10-fold decrease in monsoon precipitation, the long- 
term rate of erosion (as indicated by apatite fission-track dating) remains steady (and rapid) in 
the northern, arid part of the study area. To examine why this might happen and in particular 
how rivers with less rainfall could incise just as rapidly in the dry area as in the wet regions, we 
extracted all 3-7 km 2 catchments in the DEM, measured their channel slopes, combined them 
with a calibrated rainfall model to calculate runoff, and then calculated changes in specific 
stream power. We discovered that about half of the “erosion energy” deficit was made up by 
narrowing of channels in the arid areas. These results were reported in an attached Nature article: 
Burbank, D. W., A. E. Blythe, J. Putkonen, B. Pratt-Sitaula, E. Gabet, M. Oskin, A. Barros, and 
T. P. Ojha (2003), Decoupling of erosion and precipitation in the Himalayas, Nature, 426, 652- 
655. 


Overall, as a result of this grant, one article has been published, one is in press, another is 
in review, and a fourth is still being completed by a graduate student. 
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Abstract 

Despite the abundance in alpine terrain of glacially dissected landscapes, the magnitude 
and geometry of glacial erosion can rarely be defined. In the eastern Kyrgyz Range, a 
widespread unconformity exhumed as a geomorphic surface provides a regional datum with 
which to calibrate erosion. As tectonically driven surface uplift has progressively pushed this 
surface into the zone of ice accumulation, glacial erosion has overprinted the landscape. With as 
little as 500 m of incision into rocks underlying the unconformity, distinctive glacial valleys 
display their deepest incision adjacent to cirque headwalls. The expansion of north-facing glacial 
cirques at the expense of south-facing valleys has driven the drainage divide southwards at rates 
up to 2 to 3 times the rate of valley incision. Existing ice-flux-based glacial erosion rules 
incompletely model expansion of glacial valleys via cirque retreat into the low-gradient 
unconformity remnants. Local processes that either directly sap cirque headwalls or inhibit 
erosion down-glacier appear to control, at least initially, alpine landscape evolution. 

Introduction 

The remarkable coincidence of glacial snowlines with average height of high mountains 
suggests that glaciation may exert primary control on mountain-range elevation, irrespective of 
tectonic rates (Broecker and Denton, 1990; Brozovic et al., 1997; Montgomery et al., 2001; 
Spotila et al., 2004). Thus the topographic evolution that accompanies the transformation of 
nascent, unglaciated landscapes to highly glaciated ones is a key to unraveling the erosional 
component of a climate-erosion feedback system (Molnar and England, 1990). Comparative 
geomorphic studies support contentions that glaciation enhances valley erosion and topographic 
relief (Brocklehurst and Whipple, 2002; Kirkbride and Matthews, 1997; Montgomery, 2002; 
Small and Anderson, 1998). Observations of planimetric cirque expansion (Federici and 
Spagnolo, 2004; Olyphant, 1981) and truncated glacial valleys (Brocklehurst and Whipple, 2002) 
indicate cirque retreat could play a competing role by expanding valleys horizontally and 
removing high topography. The significance of cirque retreat in limiting mountain range height, 
however, has not been conclusively established. 

For tectonically active mountain ranges undergoing surface uplift, their progressive 
penetration into the zone of glaciation can provide a proxy for transformation of landscapes 
during the well-documented late Cenozoic cooling (Shakleton and Opdyke, 1973). Growth of the 
Tien Shan (Fig. 1, inset), when combined with its unique geologic histoiy, establishes a natural 
laboratory to investigate glacial valley erosion and cirque retreat. Basement-cored uplifts of the 
northern Tien Shan exhume resistant Paleozoic plutonic and metamorphic bedrock from beneath 
unconsolidated Cenozoic non-marine sedimentary strata. A regionally extensive, pre-Cenozoic 
low-relief unconformity separates these units (Chediya, 1986). In the arid environment of 
central Asia, significant areas of this exhumed unconformity form geomorphic surfaces that 
appear little modified by fluvial erosion, and instead faithfully delineate the geometry of range- 
scale folding (Burbank et al., 1999). When surface uplift causes ranges to intersect the local 
Pleistocene glacial equilibrium-line altitude (ELA), however, glacial erosion may cause 
significant incision of the Paleozoic rocks underlying the unconformity. Reconstruction of this 


datum over glaciated topography thus enables reliable estimation of the magnitude of glacial 
bedrock erosion. 

Late Cenozoic uplift of the Kyrgyz Range provides a structurally well-controlled 
experiment from which to measure sequential glacial landscape evolution. Uplift occurs via 
hangingwall folding and slip on a reverse fault system that emplaces resistant Paleozoic 
quartzite, metaconglomerate, granitic, and metavolcanic rocks northward over Miocene through 
Quaternary foreland-basin strata. Bedrock thermochronology supports erosion rates of ~1 
km/Myr from 3 Ma to present in the central, highest part of the range (Bullen et al., 2001). 
Detrital apatite fission-track ages from foreland-basin strata and modem rivers indicate 
progressively less exhumation in the eastern Kyrgyz Range (Bullen et al., 2003), consistent with 
preservation of exhumed, tilted remnants of the unconformity on over _ of the easternmost 50 
km of its south-facing slope (Fig. 1). These observations support an eastward propagation or an 
eastward-declining rate of coeval rock and surface uplift that has progressively elevated the 
range crest into the zone of glaciation above the local Pleistocene ELA. 

Glacial Erosion of the eastern Kyrgyz Range 

Exhumation of the pre-Cenozoic unconformity as a geomorphic surface underpins the 
quantification of erosion. Swath topographic profiles (Fig. 2) reveal minimal topographic relief 
development on its unglaciated south-facing slope. Maximum elevations of these areas are 
defined by low relief, upland surfaces with concordant elevation and southward tilt preserved 
between larger, incised drainages that descend from the range crest (Fig. 3). These surfaces 
merge with outcrops of the unconformity beneath Cenozoic strata at the base of the range. A 
fourth-order polynomial surface with curvature of <l°/km fit to the distribution and tilt of 
mapped remnants defines a continuous envelope over the south-facing slope of the range (Fig. 
4). This envelope captures the mapped unconformity within ±100 m (Ict), over >1500 m of 
structural relief. Tilting of the restored unconformity surface increases both with elevation and 
from east to west, consistent with a gradient of slip over a curviplanar reverse fault at depth. 

Comparison of three adjacent catchments incised into the south-facing slope of the 
Kyrgyz Range reveals details of the form of initial glacial erosion and the continuity of 
interfluvial geomorphic surfaces (Fig. 5). All three streams descend at steep (>8°) slopes incised 
<600 m into the tilted unconformity preserved at concordant undissected interfluves. The fluvial 
part of each valley is analogous to a half-pipe, with uniform valley width and fairly uniform 
incision. In contrast, the glaciated parts of these catchments display significant widening and 
deepening with the greatest incision at their glaciated headwaters. Adjacent glaciated catchments 
show a similar pattern of maximum incision at their headwaters (Fig. 4). 

North-facing glaciated catchments of the eastern Kyrgyz Range consist of high-elevation, 
gently sloped U-shaped valleys linked to steep, V-shaped canyons that descend to the range- 
bounding reverse fault (Fig. 1). Absence of preserved remnants of the unconformity near the 
fault limits precise estimation of incision of these catchments. To conservatively evaluate glacial 
incision by north-facing catchments, the slopes of north-facing glaciated valleys are examined 
with respect to slope of the adjacent, south-dipping unconformity surface preserved at the range 
crest (Fig. 4). The difference of these slopes defines a ratio of valley length to valley incision 
(0.3-0.4) that is significantly less than the average slope (0.8 ±0.1 (lo), n=66) of cirque 
headwalls along the range crest. 

The presence and extent of glaciation in north-facing catchments corresponds to changes 
of trend and elevation of the range crest. Outside of the zone of glaciation, the range crest tracks 
~5 km south of the trace of the range-bounding reverse fault, the range shape reflects growth of 
the underlying geologic structure (e.g. Fig. 2a), and peak elevations climb steadily from 2400 m 
near the Chu River in the east to 3400 m at the first glaciated catchment. Within the zone of 


glaciation, peak elevations stabilize at -4000 m elevation and the position of the range crest 
shifts systematically southward relative to the range front in proportion to the degree of Late 
Pleistocene ice cover (e.g. Figs. 2b through 2d). 

Discussion 

The gradient in rock uplift of the eastern Kyrgyz Range permits along-strike space-for- 
time substitution of observations of incision of the range by glaciers. Because the range is rising 
above an active reverse fault at approximately 1 km/Ma, glaciation of the eastern range tip is 
likely to have occurred sequentially during the Pleistocene. Although each valley has a unique 
incision history, spatial trends of landscape evolution that are shared among several adjacent 
valleys may represent temporal trends within an individual valley. Much of the transformation of 
the Kyrgyz Range into its glacially sculpted form involves competition between adjacent 
catchments that drives migration of drainage divides over time. Exhumation of the pre-Cenozoic 
unconformity as a geomorphic surface datum affords a unique opportunity to quantify the degree 
of glacial incision that accompanies progressive changes in valley form. 

Rapid formation of gently sloped upland U-shaped valleys accompanies initial glacial 
erosion of the eastern Kyrgyz Range. Glaciers on south-facing slopes accomplish this by 
widening and eroding their valleys deeper than adjacent fluvial valleys, evidenced by incision 
and valley-width maxima near the range crest (Fig. 5). Approximately 0.5 km of incision into the 
unconformity surface forms these characteristically glaciated valleys, consistent with low erosion 
depths required to obtain stable glacial valley cross-sections in numerical simulations (Harbor, 
1992). 

Systematic, relative southward migration of the drainage divide of the eastern Kyrgyz 
Range shows that cirque retreat also accompanies initial glacial erosion. Expansion of north- 
facing valleys via glacial cirque retreat provides a consistent explanation for these spatial trends 
of landscape evolution with progressive rock uplift above the local ELA (Fig. 6): (1) range crest 
elevation stabilizes at -4000 m with the onset of erosion by glaciers; (2) further uplift and tilting 
of the range is accompanied by southward relative migration of the drainage divide; and (3) the 
pre-Cenozoic unconformity is preserved as a geomorphic surface at the range crest. Peak heights 
remain constant along trend and the range crest preserves remnants of the exhumed 
unconformity because the divide has shifted southward with progressive uplift and glacial 
erosion (Fig. 2). In effect, during initial glacial erosion, the range crest remains more or less 
pinned to the intersection of the 4000-m contour line with the tilted unconformity. This 
explanation requires that a competitive advantage for cirque erosion exists for north-facing 
catchments over their south-facing counterparts and that glacial erosion of the Kyrgyz Range 
occurs, at least initially, by a combination of cirque retreat and glacial valley incision. 

Alternative origins for southward migration of the divide include differential fluvial 
incision, exposure of lithologies that resist glacial erosion, or formation in place of the present 
divide by range-scale folding. Range-scale folding may be ruled out as an explanation because 
peaks and ridges north of the present divide are locally at the same elevation or higher than the 
exhumed unconformity preserved at the divide (Figs. 1, 2c, and 2d). This geometry requires that 
the structural axis of the Kyrgyz Range trends north of the present drainage divide in the zone of 
glacial erosion. Metasedimentary rocks that underlie the glaciated area strike transverse to 
valleys and dip steeply, inconsistent with lithologic control of valley floor elevation. It is also 
unlikely that the present low-gradient upland valleys formed entirely by glacial incision of 
steeper fluvial valleys (c.f. Brocklehurst and Whipple, 2002) because the valley heads terminate 
at uneroded remnants of the unconformity. Only headward erosion of low-gradient glaciated 
valleys by cirque retreat could maintain this geometry at the range crest (Fig. 6). To estimate the 



total amount of cirque retreat and valley incision requires additional knowledge of the divide 
position that developed by range-scale folding and fluvial incision prior to glaciation. A tentative 
reconstruction of the original divide position across the easternmost four glaciated catchments 
indicates between 0.9 and 4.4 km of cirque retreat with only 0.4 to 1.3 km of valley incision (Fig. 
4). Further extrapolation of the pre-glaciation divide suggests that >10 km cirque retreat 
dominates erosion of north-facing valleys. 

Evidence for substantial headward expansion of catchments via cirque retreat in the 
Kyrgyz Range, the Sierra Nevada (Brocklehurst and Whipple, 2002), and potentially elsewhere 
(Olyphant, 1981) compels a closer examination of the resulting valley form. The immediate 
cause of cirque retreat is the collapse of steep headwalls at the expense of adjacent low-relief 
uplands. Processes that sap the cirque walls at the head of the glacier ultimately determine the 
retreat rate. Because the ratio of valley length to valley incision (0.3-0.4) in the headwaters of 
north-facing catchments is 2 to 3 times less than the slope of adjacent cirque headwalls (0.8 ± 
0.1), sapping by valley incision alone is unlikely to cause the observed retreat. Rather, erosion 
localized at the base of cirque headwalls is required to propagate cirque retreat across the range 
at a rate 2 to 3 times greater than the glacial valley incision rate. Incision maxima developed at 
the heads of south-facing glaciated valleys (Fig. 4) are consistent with a localized sapping 
process at cirque headwalls. 

Because cirque headwalls form within the glacial accumulation zone, modeling cirque 
erosion via ice-flux based rules (Hallet, 1979, 1996) predicts in greater downstream valley 
incision than erosion of the cirque floor: a prediction incompatible with observations from the 
Kyrgyz Range. This may be reconciled if subglacial water pressure fluctuations within cirques 
locally enhance erosion (Hooke, 1991), or if development of glacial bedforms, such as terminal 
moraines and overdeepenings, restricts pressure fluctuations and erosion down-glacier (Alley et 
al., 2003). It is also possible that the integrated residence time of vestigial cirque glaciers under 
average Quaternary conditions (Porter, 1989) has been sufficient to tip the balance of erosion in 
favor of cirque incision and consequent headward retreat. Such erosion is consistent with 
preferred erosion of north-facing catchments, where shading of glacier surfaces lowers the local 
ELA by ~200 m. 

Conclusion 

Analysis of the configuration and gradient of glacial valleys in the eastern Kyrgyz Range 
indicates that cirque retreat dominates initial glacial erosion. This erosion is quantified in the 
context of a progressively tilted pre-Cenozoic unconformity that is exhumed as a geomorphic 
surface. This cirque retreat may be incompatible with existing ice-flux-based glacial erosion 
rules, and highlights a need to understand erosional processes localized within cirques in order to 
predict alpine landscape evolution. Cirque retreat can effectively bevel across an elevated alpine 
plateau and keep pace with moderate rock uplift rates, as documented here for the Kyrgyz Range 
where retreat rates are at least 2 to 3 times greater than the accompanying valley incision rate. 
Although cirque retreat may not be as significant in all alpine settings (White, 1970), results 
from the Kyrgyz Range support the contention that glacial valley incision combined with cirque 
retreat into interfluves may be a viable mechanism whereby glaciers simultaneously erode at 
high rates and limit alpine relief. 


Figures 



Exhumed 
Unconformity jt 


Figure 1. Exhumed pre-Cenozoic unconformity and maximum glacial extent in the eastern 
Kyrgyz Range. Present-day drainage divide shown as line of peaks (elevations in meters) 
connected by heavy gray solid line. A-D are center lines of swath profiles in Fig. 2. Inset: shaded 
relief map of central Asia with Kyrgyz Range crest as hatched line. 
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Figure 2. Topographic swath profiles across the easternmost Kyrgyz Range. Each profile 
samples topography from a 5-km-wide swath centered on profile lines shown on Figure 1 . Light 
gray band depicts range of minimum and maximum topography, with mean values shown as 
thicker central line. Horizontal black arrow depicts extent of geomorphic surfaces formed by 
exhumed remnants of the pre-Cenozoic unconformity'. ELA is last glacial maximum equilibrium- 
line altitude. White arrow is extent of last maximum glaciation. Glaciated areas correspond to 
areas of higher relief and erosional removal of the pre-Cenozoic unconformity at the range crest. 
Extent of glaciation is inversely proportional to northward position of drainage divide, illustrated 
by vertical dashed lines. Highest elevation of the range crest shown by downward arrows. 



Figure 3. A. ASTER scene of geomorphic surfaces formed by exhumed pre-Cenozoic 
unconformity on south slope of the Kyrgyz Range. Surfaces are visible as medium albedo with 
intricate dendritic stream network. B. Interpreted slope map of area shown in A. Steeper slopes 
(darker gray) correspond to valley sides. White solid lines outline surface remnants mapped from 
analysis of the SRTM topography, ASTER, and Corona imagery (not shown). Elevation contours 
at 100-m intervals (black lines) show minimal local relief formed by incision of geomorphic 
surfaces. White dashed line shows folded Paleozoic contact (truncated by the unconformity) 
between metavolcanic rocks (mv.) and quartzite (q). 











Figure 4. Contours of reconstructed pre-Cenozoic unconformity elevation (solid lines) and slope 
(dashed lines) overlain on a shaded relief image of the Eastern Kyrgyz Range derived from 
SRTM topography. AUnconformity depicts elevations within 100 m and 200 m from modeled 
unconformity and compares well with its mapped exhumed remnants on Fig. 1. Dotted lines 
show gradients of north-facing glacial valley floors. Pre-glaciation divide extrapolated along 
front of prominent faceted spur ridges that separate glaciated valleys. 


Figure 5. Valley profiles, incision beneath the reconstructed unconformity, and catchment width 
from three adjacent south-facing streams with varying amounts of glaciation. See Fig. 4 for 
catchment locations. In the absence of glaciation, incised valleys (SI) are shaped like half-pipes 
with uniform width. Glacial erosion corresponds to incision maxima at cirque floors and 
widening of catchments in the zone of ice accumulation, supporting valley expansion via cirque 
retreat. 
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Figure 6. Landscape evolution of the eastern Kyrgyz Range with progressive rock uplift. 

Compare to Figs. 2a through 2c. Light and dark gray represent ridgeline and valley elevation, 
respectively. Dashed line shows extrapolation of folded unconformity surface over north slope of 
range. 1. Range form prior to glacial erosion with divide controlled by fluvial erosion and range- 
scale structural fold axis. 2. Initial glacial erosion of north-facing valley causes southward divide 
migration relative to fold axis. Unconformity surface is preserved on south-facing slope at range 
crest. 3. Further range uplift accompanied by additional cirque retreat and expansion of north- 
facing valley. Note peak of northern ridgeline is higher than unconformity preserved at the 
drainage divide. 
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Effects of Bedrock Landsliding on Cosmogenically Determined Erosion Rates 
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The successful quantification of long-term erosion rates underpins our understanding of landscape formation, the topographic 
evolution of mountain ranges, and the mass balance within active orogens. The measurement of in situ-produced cosmogenic 
radionuclides (CRNs) in fluvial and alluvial sediments is perhaps the method with the greatest ability to provide such long- 
term erosion rates. In active orogens, however, deep-seated bedrock landsliding is an important erosional process, the effect 
of which on CRN-derived erosion rates is largely unquantified. We present a numerical simulation of cosmogenic nuclide 
production and distribution in landslide-dominated catchments to address the effect of bedrock landsliding on cosmogenic 
erosion rates in actively eroding landscapes. Results of the simulation indicate that the temporal stability of erosion rates 
determined from CRN concentrations in sediment decreases with increased ratios of landsliding to sediment detachment rates 
within a given catchment area, and that larger catchment areas must be sampled with increased frequency of landsliding in order 
to accurately evaluate long-term erosion rates. In addition, results of this simulation suggest that sediment sampling for CRNs is 
the appropriate method for determining long-term erosion rates in regions dominated by mass-wasting processes, while bedrock 
surface sampling for CRNs is generally an ineffective means of determining long-term erosion rates. Response times of CRN 
concentrations to changes in erosion rate indicate that climatically driven cycles of erosion may be detected relatively quickly 
after such changes occur, but that complete equilibration of CRN concentrations to new erosional conditions may take tens of 
thousands of years Simulation results of CRN erosion rates are compared with a new, rich dataset of CRN concentrations from 
the Nepalese Himalaya, supporting conclusions drawn from the simulation. 


1. Introduction 

The successful quantification of long-term ero- 
sion rates underpins our understanding of landscape 
formation, the topographic evolution of mountain 
ranges, and the mass balance within active orogens 
[1]. Observed changes in long-term erosion rates are 
often considered proxies for changes in the climatic or 
tectonic boundary conditions that control landscape 
evolution. As the importance of long-term erosion 
rates to unraveling significant geomorphic and tec- 
tonic problems has become clear, a widely applica- 
ble means of measuring long-term surface denudation 
rates has been sought. The measurement of in situ- 
produced cosmogenic radionuclides (CRNs) in fluvial 
and alluvial sediments has been shown to yield spa- 
tially averaged erosion rates, and has become perhaps 
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the method with the greatest applicability in measur- 
ing erosion rates over 1 0 3 — 1 0 5 years and across a 
wide variety of landscapes and erosional processes 
|2— 6]. Most existing studies that utilize CRN -derived 
erosion rates have focused on regions with fairly 
spatially and temporally homogeneous erosion rates. 
Here we present a simulation to explore the effects of 
bedrock landsliding on cosmogenic erosion rates, and 
potential for exploiting cosmogenic nuclides to mea- 
sure erosion rates in rapidly eroding, active orogens. 

Large bedrock landslides incise to depths greater 
than one or more attenuation lengths of cosmic rays, 
thus mobilizing sediments with little or no cosmo- 
genic nuclide abundance (Fig. 1). For example, based 
on empirical relationships of landslide depth to area 
17], a landslide with a radius of just 10 meters will in- 
cise to ~ 100 cm, below one attenuation depth for 
spallogenic nuclide production. An extraordinarily 
large slide may incise to a maximum depth of sev- 
eral tens of meters, below one attenuation depth for 
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Figure 1. Cartoon showing the effects of landslides on surface cosmogenic nuclide concentrations. A) Homogeneous CRN surface concen- 
tration and sediment volume and concentration during steady-state erosion. B) Heterogeneous surface CRN concentration and heterogeneous 
sediment volume and CRN concentration during landslide-dominated erosion. 


muogenic nuclide production. In catchments where 
deep landslides have recently occurred, the addition 
of nuclide-poor landslide detritus to the fluvial system 
will dilute the CRN concentration in the fluvial sed- 
iment, yielding apparently higher erosion rates. For 
example, samples from low-order fluvial catchments 
in the Nepalese Himalaya yield ‘erosion rates’ de- 
termined from cosmogenic nuclides that range from 
0.01 to 0.5 mm/yr, even along the same ridge crest 
(Heimsath, in prep.). Similar results are reported from 
the San Bernardino Mountains in southern California, 
where ‘erosion rates’ along a landslide-dominated es- 
carpment vary from < 0.3 to > 2.7 mm/yr, while rates 
from an adjacent region dominated by sediment de- 
tachment are on the order of hundredths of mm/yr 
18]. Whereas such variability may be expected for 
data sets that are focused on small, steep catchments 
in actively deforming mountain belts, it does not ex- 
plain how in situ-produced CRNs may be utilized to 
understand the rates of erosion in regions affected by 
both sediment detachment and landsliding processes. 

Here we address the effects of bedrock landslid- 
ing on CRN-derived erosion rates from bedrock and 


fluvial samples in an attempt to answer several basic 
questions: As the rate of landsliding increases, how 
are CRN concentrations in fluvial sediment affected? 
Can reliable CRN erosion rates be derived from flu- 
vial sediment or bedrock samples when landsliding 
is the dominant erosional process? Over what spatial 
scales do fluvial systems integrate the effects of land- 
sliding? What timescales are required for CRN con- 
centrations to respond to changes in erosion rates? We 
present the results of a numerical simulation of cos- 
mogenic nuclide production and erosional removal in 
landslide-dominated catchments to assess these ques- 
tions. A series of simulations with varying sediment 
detachment and landslide erosion rates are used to 
create statistical populations of CRN-derived erosion 
rates for both ‘sediment’ and ‘bedrock’ samples for a 
theoretical landscape. Finally, simulated distributions 
of CRN-derived ages are compared to a new', rich data 
set of CRN-derived erosion rates from the Nepalese 
Himalaya. 


Effects of Bedrock Landsliding on Cosmogenic Erosion Rates 


3 


2. Numerical Simulation 

The numerical simulation is based on actual digital 
elevation data, and simulates the production of cos- 
mogenic nuclides at each model cell, the removal of 
material through sediment detachment, the removal 
of material by landsliding, and the radioactive de- 
cay of cosmogenic nuclides. A Geographic Informa- 
tion System (GIS; in this case Arclnfo) is used for 
the backbone of the simluation. Model initialization, 
data assimilation, and data output are all controlled 
through the GIS using Arc Macro Language (AML). 
Computationally intensive portions of the model are 
passed from the GIS to customized Perl modules for 
computational efficiency. The three main functional- 
ities of the simulation are described in greater detail 
below. 


2.1. Cosmogenic Nuclide Production 

Prior to running the landslide simulation model, 
cosmogenic nuclide production rates must be calcu- 
lated for each cell in the model. Calculation of these 
rates begins with a geo-referenced digital elevation 
model (DEM) of the study area of interest. For each 
cell in the digital elevation model, cosmogenic pro- 
duction scaling factors are calculated within the GIS 
based on cell altitude and latitude following [9,10]. 
Further corrections to cosmogenic production are ap- 
plied by calculating the topographic shielding at each 
point in the DEM. For each cell in the DEM, the ver- 
tical angle to to every other cell is calculated. These 
values are binned into 5° radial bins, and the maxi- 
mum vertical angle in each bin is used to approximate 
the horizon angle for that bin. The topographic shield- 
ing factor for each bin is derived from the horizon 
angle using a published methodology [11]. The al- 
titude, latitude, and topographic shielding factors are 
combined within the GIS system to produce an output 
array of cosmogenic production scaling factors. This 
array is then multiplied by the high-latitude, sea-level 
production rate of the cosmogenic nuclide of interest 
to create an array of cosmogenic production rates. In 
this case, we have chosen to model J0 Be and selected 
a production rate of 5.3 atoms/g/yr [12]. This array of 
10 Be production rates is preserved for use through the 
rest of the model run. 


22. Model Initialization 

Model initialization consists of two separate ac- 
tions: i) preparing the model for data gathering and 
assimilation and ii) calculating an initial surface cos- 
mogenic nuclide concentration to start the landslide 
model. The first of these two tasks is the most time- 
intensive and must be performed separately for each 
DEM on which this model is run. Using the hydro- 
logic functions available in a GIS, watersheds within 
the model area are delineated. A variety of first- 
through highest-order watersheds are selected and 
saved for later use in data analysis. Additionally, 100 
random points are generated across the model space 
at which to track cosmogenic nuclide concentrations 
in bedrock and erosional removal of material. After 
the data gathering initialization steps are complete, an 
input cosmogenic surface nuclide concentration array 
must be calculated for input to the landslide portion of 
the model. We have chosen to input a surface concen- 
tration grid that represents the steady-state concentra- 
tion of 10 Be at the sediment detachment rate specified 
for the model run. The sediment detachment rate is 
limited by the rate at which rock can be converted to 
soil or regolith, and is taken to be less than 0.3 mm/yr 
[13]. The surface concentration ( N ) is calculated for 
each point in the mode! following this equation: 


N 10 Be = 


P 10 Be 

(Aio Be -f E/A) 


(1 _ e ~(^Be + E / A '> t ) (]) 


Where jV 10 Be and P 10 Be are the concentration 
(atoms/g) and production rate (atoms/g/yr) of 1() Be, 
respectively; Aio Bc is the decay constant of 10 Be 
(yr _1 ); E is the erosion rate (here equivalent to the 
sediment detachment rate), in g/cm 2 /yr; A is the neu- 
tron attenuation length in rock, in g/cm 2 ; and t is time 
(yr). This initial concentration array is saved for input 
into future model runs. 


2 _3. Landslide Simulation 

The landslide simulation portion of the model takes 
as input the surface cosmogenic nuclide concentra- 
tion array calculated in the previous step, and the cos- 
mogenic production rate array calculated in the first 
model step. This portion of the model is iterated. 
At the end of each iteration, two arrays are output, a 
depth array which contains the sum of all sediment 
removed by erosion, due to both sediment detach- 
ment and landsliding, and a surface concentration ar- 
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ray that contains the surface cosmogenic nuclide con- 
centration at each model cell after erosional removal 
of material, cosmogenic ingrowth and radioactive de- 
cay. The individual steps are detailed below. 

23.1. Cosmogenic ingrowth and decay 

The first step in the iteration of the landslide model 
is to add cosmogenic nuclide ingrowth for one model 
time step (t) to the input surface cosmogenic nuclide 
concentration array. Given, for each cell in the model, 
an initial surface cosmogenic nuclide concentration, 
./V t 10 Be, and a surface cosmogenic nuclide production 
rate, P 10 Be, a resultant surface cosmogenic nuclide 
concentration, AfJ°Be, produced by cosmic ray bom- 
bardment and removal by radioactive decay can be 
calculated by: 

AT r 10 Be = (iV/ 0 Be + P 10 Be • t)e~ Xl0 B «* (2) 

The new surface concentration, 7V r 10 Be, replaces 
the initial value in the surface concentration array. 

2.3.2. Sediment detachment 

For each model run, a sediment detachment rate, 
E s is specified. This rate represents the spatially ho- 
mogeneous erosional removal of material from the 
land surface, limited by weathering and soil produc- 
tion processes (e.g. grussification in granitic ter- 
ranes). For each model time step (t), a depth equal 
to E s x t is removed from the landscape and added to 
the depth grid. This material is always removed from 
the upper surface of the topography, and therefore has 
the highest concentration of CRNs. 

2.3.3. Landslides 

After removal of material by sediment detach- 
ment, the model is populated with landslides. Land- 
slides are assumed to obey a power-law frequency- 
magnitude relationship [7,14,151. Based on this as- 
sumption, populations of landslides in the model can 
be derived from four parameters: /3, the power-law 
exponent for landslide frequency-magnitude relation- 
ship, A m i n , the minimum landslide area considered 
in the model, A max , the maximum landslide area con- 
sidered in the model, and Ei s , the average rate of ero- 
sion by landsliding over the model area. Although 
short-term landslide erosion rates fluctuate due to the 
episodicity of landslides, the average rate of ero- 
sion by landsliding is produced through the power- 
law frequency-magnitude relationship over time (Fig. 



Figure 2. Comparison of output model landslide erosion rates with 
prescribed landslide erosion rates, averaged over 10,000 year time 
steps. Mean erosion rates for the 1 00,000 year model run are shown 
at the right of the figure. Short term variations in landslide erosion 
rate are due to random distribution of slide size; long term landslide 
erosion rates return input erosion rates ±3%. 


2). Complete derivations of landslide frequency- 
magnitude distributions in the model space are given 
in Appendix A. For each model timestep, a land- 
slide distribution is generated, and the landslides are 
randomly distributed over the model space. The to- 
tal amount of material removed by landsliding from 
each model cell during the timestep is then calculated. 
The total depth of material removed by landsliding is 
added to the depth of material removed by sediment 
detachment (see Appendix A). 

2.3.4. Surface concentration 

Once the depth array is tabulated, the surface cos- 
mogenic nuclide concentration is recalculated to re- 
flect the depth of material removed from each model 
cell. Assuming, for simplification, a constant rock 
density ( p ) in the model of 2.65 g/cm 8 , and a cos- 
mogenic ray attenuation length in this rock (A) of 150 
cm 2 /g, a final surface cosmogenic nuclide concentra- 
tion at the end of the timestep, Nj° Be. can be cal- 
culated for each cell based on the cosmogenic con- 
centration following ingrowth, A r r 10 Be. above, and the 
depth of material. D. removed from the cell during the 
timestep: 

A r y°Be = A r r 10 Be ■ e~ Dip/A '> (3) 

The value Nj°Be is stored for each cell in the surface 
concentration array. At this point, this array is saved 
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for data extraction, as described in the next section, 
and then re-opened as the initial surface concentration 
grid for the next iteration step of the landslide portion 
of the model. 


2.4. Data Extraction 

At a user-specified sampling interval, the depth ar- 
ray and surface concentration array are passed back 
to the GIS, and data are extracted from each array to 
simulate two potential cosmogenic nuclide sampling 
methods: surface exposure age dating and stream sed- 
iment sampling. 

2.4.1. Surface exposure age dating 

To simulate the effects of landsliding on surface 
age exposure dating, at each sampling interval, the 
depth of material removed from, and the final sur- 
face concentration at, each of 100 randomly gener- 
ated points across the model are recorded. From these 
data, a volumetric erosion rate ( E v ) for this sampling 
interval is calculated by dividing the depth of material 
removed ( D ) by the model time step (f; recall that D 
includes both continual removal of rock by sediment 
detachment and episodic landsliding). Additionally, a 
cosmogenic erosion rate ( E c ) is calculated, for each 
point using the standard assumption of steady-state 
erosion [16]: 


E c 


A ( P 10 Be 
p \N { W Be Al ° 


(4) 


Both volumetric and cosmogenic surface-exposure- 
age erosion rates, along with the depth and surface 
concentration information, are stored in a database 
file. Data can be extracted after the completion of 
the model run to analyze the variation in volumetric 
and cosmogenic erosion rates at a given point through 
time, or data from all points can be assimilated and a 
probability density function (PDF) for either erosion 
rate can be calculated. 


2.4.2. Stream sediment sampling 

In addition to surface age exposure dating, cosmo- 
genic nuclide concentrations in stream sediments can 
be used to estimate average upstream erosion rates. 
To test the effects of landsliding on such estimates of 
erosion rates, the model also extracts data by water- 
sheds. At each sampling interval, the following pro- 
cedures are carried out for each watershed of inter- 
est. A volumetric erosion rate (E v ) is calculated by 


summing the total depth of material removed from the 
watershed and dividing it by the watershed area mul- 
tiplied by the model time step. Second, a cosmogenic 
erosion rate ( E c ) for the watershed is calculated. To 
calculate this erosion rate, first the concentration of 
the cosmogenic nuclide of interest in the eroded sedi- 
ment must be determined. Using the depth array and 
surface concentration array, the average nuclide con- 
centration in the eroded material at each model cell, 
A r l, 0 Be, can be calculated as follows: 

AT^Be = AT/ 10 Be • A (l - e - D(p/A) ) (5) 

The concentration of cosmogenic nuclides in the sed- 
iment removed from the watershed, AT| 0 Be, then, is: 

^ 0 Be = ^A“Be/^P (6) 

The erosion rate derived from the cosmogenic nuclide 
concentration in the material removed from the water- 
shed is then calculated as in Eqn. 4, where P 10 Be 
would represent the watershed-averaged production 
rate (derived from the GIS-based production calcula- 
tions), and A^°Be would be replaced with Ng°Be. 
Both volumetric and cosmogenic erosion rates are 
recorded to database files, and saved for later anal- 
ysis. 

2.43. Simplifications and assumptions 

It should be noted that the model outlined above 
makes several simplifying assumptions regarding sed- 
iment production and transport. With regards to sedi- 
ment detachment and transport, this is assumed to oc- 
cur at a homogeneous rate across the entire model, 
although in reality the rate at which these processes 
act on the scale of our model are controlled by lo- 
cal slope and lithology. Landslides are also randomly 
placed across the landscape, with no considerations 
for hillslope and aspect. Second, the model contains 
no provision for sediment storage. All material de- 
rived from landsliding is assumed to pass through 
the model within the timestep in which the landslide 
occurred (100 years for the models discussed here). 
Third, the model is intended to produce a population 
of cosmogenic and volumetric erosion rates for statis- 
tical analysis. Erosion and landslides occur through 
‘time’ to produce a variety of surface CRN concen- 
trations that could potentially be sampled; however, 
this model is not a landscape evolution model. Dur- 
ing the course of the model run, the model landscape 
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Table 1 

Statistics for basins tracked in the model. 


Basin Order 

Number of 
Basins 

Mean Area 
(km 2 ) 

First 

30 

0.1 

Second 

6 

0.6 

Third 

5 

2.3 

Fourth 

4 

8.4 

Fifth 

1 

69.0 


surface does not evolve, and shielding effects or ab- 
solute elevation changes that in reality would alter 
CRN production rates are not considered. Finally, the 
role of muogenic production is not explicitly consid- 
ered in this simulation. The greater attenuation length 
of muons would likely have the effect of moderating 
to some extent the surface nuclide concentration in 
the wake of small- to moderate-sized landslides, but 
large landslides, which carry the majority of the ero- 
sional load, will incise more deeply than the attenua- 
tion depth of muons, and thus would have a negligible 
effect on the model. 

3. Model Results 

To examine the model, we simulated the effects 
of landsliding on cosmogenic nuclide equilibrium, 
and associated CRN-derived and volumetrically cal- 
culated erosion rates for the San Antonio Creek catch- 
ment, located in the eastern San Gabriel Mountains of 
southern California (Fig. 3). It is a small (~ 70 km 2 ), 
mountainous catchment, selected in part because of 
the availability of high-quality digital elevation model 
(DEM) data over the region (30-meter resolution), 
and in part because a significant amount of work ex- 
ists describing the geomorphology and neotectonics 
of the region. Low -temperature thermochronologic 
data indicate that the eastern San Gabriel Mountains 
are being exhumed at a rate of ~ 0.3 — 1 mm/yr f 17- 
19], while geomorphic and geologic studies indicate 
that landsliding is a prevalent mechanism of erosion 
in this watershed [20,15]. In fact, both a landslide 
frequency-magnitude exponential scaling factor (/?) 
and a long-term average erosion rate have been de- 
termined for this region [15]. Interestingly, in the San 
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Figure 3. Shaded relief map of the San Antonio Creek watershed. 
San Gabriel Mountains. California. The full watershed is outlined 
with a heavy black line, while representative smaller order water- 
sheds are indicated with a lighter black line. Random ‘bedrock’ 
sampling points (discussed in the text) are shown as black dots. 


Gabriel Mountains, J3 = 1.1, similar to results from 
the Southern Alps of New Zealand [7]. When 0 is less 
than 1.5, large, but infrequent landslides dominate the 
overall sediment flux from a catchment. As such, this 
watershed potentially provides a natural laboratory to 
study the effects of landsliding on CRN-derived ero- 
sion rates, and for comparison with and calibration of 
the numerical model. 

Within the San Gabriel Mountains, San Antonio 
Creek is a fifth-order stream. The San Antonio Creek 
watershed was divided into 46 sub-basins that were 
tracked as part of this simulation, from which mod- 
eled sediment-derived CRN erosion rates were calcu- 
lated (examples shown in Fig. 3; basin statistics are 
listed in Table 1). In addition to the sub-basins, 100 
points were randomly distributed across the model 
space to serve as simulated bedrock CRN sampling 
localities (Fig. 3). For San Antonio Creek, sedi- 
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Figure 4. Example of a box -and- whisker plot to graphically display 
the statistical distribution of a non-Gaussian data set. A probability 
density function (PDF) and boxplot are shown for actual model 
results of cosmogenic erosion rates from first-order basins with a 
sediment detachment rate of 0.1 mm/yr and a total erosion rate of 
10 mm/yr. The distribution of erosion rates is plotted across the 
bottom of the graph as a series of vertical black hash marks. 


ment detachment rates of 0.01, 0.1 and 0.3 mm/yr 
were selected, spanning the range of observed sedi- 
ment detachment rates in the Transverse Ranges [8] 
to the most rapid known rates of rock-to-regolith con- 
version [13]. For each of the three sediment detach- 
ment rates, simulations of total erosion rates of 1, 
5, and 10 mm/yr were run. The rate of erosion by 
landsliding in each simulation run is the difference 
between the total erosion rate and the sediment de- 
tachment rate. Once the simulation was equilibrated 
to the imposed erosion rate, the simulation ran for 
100,000 yrs in 100 yr timesteps. Simulated cosmo- 
genic nuclide concentrations were sampled from the 
46 tracked catchments and the 100 bedrock locations 
every 1000 years, yielding 100 data points per analy- 
sis element (either catchment or bedrock location). 

Results from the simulation are summarized and 
displayed as a box-and-whisker plot. Such a plot al- 


lows a fairly straight-forward visual presentation of 
the statistical distribution of a data set, and is a valu- 
able graphical method to compare the distribution of 
two separate data sets (Fig. 4). 

3.1. Simulated sediment erosion rates 

Catchment-wide erosion rates from the 46 catch- 
ments tracked in the model allow comparison of 
the statistical distributions between CRN-derived and 
volumetrically averaged erosion rates for each of the 
nine simulations, and illustrate variations in the statis- 
tical distributions within each simulation as a function 
of catchment order (Fig. 5). 

Several conclusions can be drawn from these data. 
First, with increasing proportion of sediment detach- 
ment, the CRN-determined erosion rate at any catch- 
ment scale more closely reflects the volumetric ero- 
sion rate for any given total combined erosion rate 
(Fig. 6). This is not unexpected, as sediment de- 
tachment is modeled as a uniformly continuous pro- 
cess in the simulation, such that increased ratios of 
sediment detachment-to-landsliding will result in a 
greater contribution from a steady-state process to 
the overall erosional volume. (It is worth noting 
that under this formulation, for any given total ero- 
sion rate, an increase in sediment detachment rate re- 
sults in a decrease in the rate of erosion due to land- 
sliding. This, increasing the rate of sediment de- 
tachment CRN-derived erosion rates that are closer 
to the volumetric erosion rates, but further from the 
total erosion rate due the reduction in frequency of 
landslides; Fig. 6). Additionally, the data empha- 
size that the observed magnitude-frequency relation- 
ships of landslides skew volumetric erosion towards 
the larger, more infrequent, landslides [14,15], This 
is reflected in all nine simulations, where the median 
50% of the CRN-derived and volumetric erosion rates 
generally fall below the total imposed erosion rate in 
the simulation. The outer statistical bounds of the 
CRN-derived and volumetric erosion rates are sub- 
stantially higher than the imposed erosion rates, re- 
flecting the the infrequency, but importance, of these 
large events in controlling erosion rates in landslide 
dominated catchments. 

The damping and averaging effect of erosion rates 
derived from CRNs in sediment are also illustrated 
in the statistical spread of volumetric- versus CRN- 
derived erosion rates at all catchment scales, sediment 
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Figure 5: Box-and-whisker plots of CRN-derived erosion rates and volumetric erosion rates for nine separate simulations with three combinations of sediment detachment and total 
erosion rates. Upper solid line on each graph indicates total erosion rate ; lower solid line indicates sediment detachment rate. See text for detailed discussion. 
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detachment rates, and landslide rates, where the dis- 
tribution of CRN-derived erosion rates is tighter than 
contemporaneous volumetric erosion rates. This ef- 
fect is particularly notable at small to moderate catch- 
ment scales in regions with low rates of erosion by 
landsliding. In particular, at 1 mm/yr total erosion 
rate, the CRN-derived rates have a 50% smaller dis- 
tribution than the volumetric erosion rates. 

The effect of spatial averaging on the spread of 
volumetric- and CRN-derived erosion rates is also 
highlighted in our results (Fig. 5). The statistical 
spread of both volumetric erosion rates increases with 
increasing catchment size, presumably related to the 
increased likelihood of experiencing a large mass- 
wasting at greater catchment area, until a threshold is 
reached, at which point the catchment becomes large 
enough to adequately average large landsliding events 
and areas unaffected by mass wasting, and the spread 
of the data drops significantly. This spatial scale in 
our simulations appears to occur between fourth- and 
fifth-order catchments (a jump from ~ 8 km 2 to ~ 70 
km 2 ). The exception to this rule is at low total erosion 
rates, where landslides appear to occur infrequently 
enough at all catchment scales to never significantly 
increase the statistical distribution of erosion rates at 
third- through fifth-order catchments. 

Finally, given the damping and spatial averaging 
affects described above, it is heartening to note that 
the median 50% of observations are more or less 
consistent between volumetric and CRN-derived ero- 
sion rates (Fig. 5). At small catchment scales, 
CRN-derived rates are typically higher than volumet- 
ric rates (this is particularly clear at low total ero- 
sion rates), but converge at larger catchment scales 
(Fig. 6). The catchment size at which this con- 
vergence occurs varies, and decreases with increased 
rates of erosion. As a general rule of thumb, how- 
ever, it would appear that CRN-derived rates of ero- 
sion from sediments are statistically representative of 
volumetric rates in our numerical simulation at third- 
, or at most, fourth-order catchment scales. This 
observation indicates that CRN-derived erosion rates 
from sediments in landslide-dominated catchments 
may in fact be useful for looking at basin-wide ero- 
sion. Even modest-scale drainages, such as San An- 
tonio Creek, have several fourth-order catchments. 
Measuring CRN-derived erosion rates from each of 
these catchments should allow the identification of 





Volumetric Erosion Rate (mm/yr) 


Figure 6. Scatter plots of median cosmogenic versus median volu- 
metric erosion rates for (from top to bottom) total erosion rates of 1 . 
5, and 10 mm/yr). All plots indicate that cosmogenic erosion rates 
are higher than volumetric erosion rates, and illustrate the conver- 
gence of the two rates with increasing sediment detachment rate. 
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outlying rates that may be due to recent mass wasting 
events. Such outlier rates would approximate the re- 
cent volumetric erosion rates from the basin (Fig. 5). 
Nonetheless, even fourth-order catchments are likely 
to underpredict the overall erosion rate by 20 - 40%, 
due to the impact of rare, but unusually large, land- 
slide events. 

3.2. Simulated bedrock erosion rates 

In addition to tracking the CRN concentration at 
46 catchments within the simulation, we also calcu- 
lated the surface CRN concentration at 100 points 
randomly distributed throughout the model space, 
and analyzed such concentrations for erosion rates as 
though bedrock samples were collected from each of 
these points (eqn. 4). The statistical distributions 
of these erosion rates are derived from analysis of 
100 random points, sampled 100 times at 1,000 year 
intervals, over a 100.000 year model run, such that 
each box-and-whisker in the plot represents the sta- 
tistical distribution of some 10,000 independently de- 
termined erosion rates (Fig. 7). 

A first-order observation is that the spread of ero- 
sion rates increases with increasing rates of total ero- 
sion (and thus, with increasing rates of landsliding). 
This almost certainly reflects the increased likelihood 
of any of the randomly sampled points in the model 
to be affected by landsliding with the increasing fre- 
quency of events. A second important observation is 
that at low rates of total erosion (and landsliding), a! 
though sediment detachment encompasses, at most, 
30% of total erosion, the median CRN-derived ero- 
sion rate is almost identical to the sediment detach- 
ment rate. That is, for any given point in the land- 
scape. the likelihood of sampling a point that has been 
recently enough effected by mass wasting to alter the 
CRN concentration is virtually negligible. As the to- 
tal rate of erosion increases, the median erosion rates 
increases above the background sediment detachment 
rate, yet fall well below the total erosion rate over the 
landscape. In contrast to the distribution of sediment- 
derived CRN erosion rates (Fig. 5), the upper extent 
of the ‘whisker’ (3 inter-quartile ranges beyond the 
median) never exceeds the imposed total erosion rate. 
Although the CRN-derived erosion rates do not fol- 
low a Gaussian distribution, one may draw the anal- 
ogy that less than 99% of the ‘bedrock’ samples will 
fail to accurately reflect the total erosion rate over the 


landscape. (The highest actual percentage is 0.25% 
for a total erosion rate of 10 mm/yr, and a sediment 
detachment rate of 0.01 mm/yr, where 25 samples 
out of 10,000 fell within 2 mm/yr of the total erosion 
rate). 

These results suggest that sampling bedrock expo- 
sures in basins dominated by mass wasting may pro- 
vide an upper bound on sediment detachment rates 
across the basin, but are inadequate as a means to de- 
rive overall rates of erosion from the combined effects 
of sediment detachment and landsliding. 

33. Response of CRN-derived erosion rates to 
changes in rates of mass wasting processes 

In addition to tracking the CRN concentrations in 
sediment from catchments, and at individual points, 
the simulation also calculates the mean CRN concen- 
tration across the landscape at each point in the sim- 
ulation. Whereas the calculation of CRN and volu- 
metric erosion rates is only performed while the sim- 
ulation is in an erosional ‘steady -state’ for the im- 
posed landslide and sediment detachment rates, the 
mean CRN concentration is calculated throughout the 
model run to gather an estimate of the response time 
of CRNs in the landscape to changes in erosional 
boundary conditions (Fig. 8). 

These plots record the mean 10 Be concentration 
over the landscape, beginning at the initialization of 
the model, where the 10 Be concentration at each cell 
in the simulation is analytically solved for the im- 
posed sediment detachment rate and the scaled pro- 
duction factor at each cell. Subsequently, landslides 
are populated across the model landscape, and the 
mean concentration of 10 Be begins to decrease at a 
rate controlled by the sediment detachment rate and 
the total rate of landsliding imposed on the model. 
After a period when mean 10 Be concentration steady- 
state is achieved during sediment detachment and 
landsliding, the landslides are eliminated from the 
simulation, and the mean 10 Be concentration in the 
landscape increases (Fig. 8A). 

Tracking of specific points emphasizes the effects 
of changes in sediment detachment and landsliding 
rates on the response time of mean 10 Be concentra- 
tion (Fig. 8B). The response time of the landscape 
to achieve a new equilibrium mean 10 Be concentra- 
tion is a function of both the sediment detachment 
rate and the landsliding rate. For a given total ero- 
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Sediment detachment rate (mm/yr) 


Figure 7. Distribution of cosmogenically determined erosion rates (shaded boxes) from individual sample points for each model run. Total 
erosion rates for each set of runs are indicated by solid lines; sediment detachment rates are drawn across the lower portion of the graph 
in dashed lines. As landslide erosion rates increase, the spread and median of cosmogenically derived erosion rates also increase. Median 
cosmogenically determined erosion rates at these points, however, remain well below total erosion rates for each model run. 


sion or sediment detachment rate, the time for 10 Be 
to reach a new equilibrium decreases with an increase 
in the other rate. Increasing the sediment detachment 
rate by a factor of ten decreases the time necessary 
for 10 Be to reach a new equilibrium concentration 
(squares on Fig. 8B); likewise, increasing the rate of 
landsliding by a factor of 5 causes also yields a de- 
crease in the response time (circles on Fig. 8B). On 
the other hand, the response time for re-equilibration 
of mean 10 Be concentration at the cessation of land- 
sliding is a function solely of sediment detachment 
rate. The time necessary for a landscape to recover 
95% of its pre-landslide mean 10 Be concentration is 
as little as 2,000 yrs, at sediment detachment rates of 
0.3 mm/yr, and as great as 55,000 to 75,000 years at 
sediment detachment rates of 0.01 mm/yr (triangles 
on Fig. 8B). 

These results indicate that in rapidly eroding land- 
scapes, changes in landsliding rates over relatively 
short time periods (say a few thousand years or less) 
may generate transients in CRN concentrations last- 
ing tens of thousands of years. If landsliding is a 
more prevalent mechanism of erosion under certain 
climatic regimes (e.g. during interglacial periods; 
[21,?]), then the response time of mean 10 Be concen- 
tration to changes in landsliding rate must be consid- 


ered when the frequency of these rate changes is of 
order the 10 Be response time of tens of thousands of 
years. 

3.4. Comparison of simulation results of CRN ero- 
sion rates to CRN erosion rates from the 
Nepalese Himalaya 

The Khudi River in Nepal, a tributary to the 
Marsyandi, has been the focus of a multi-disciplinary 
study of geomorphic and geodynamic coupling in the 
Flimalaya. As part of this study, erosion rates over 
this catchment have been assessed using a variety of 
techniques, including long-term erosion rates from 
low-temperature thermochronometers (> 2-5 mm/yr; 
1 22]), and present-day erosion rates from stream- 
sediment suspended-load determinations (~ 3 mm/yr; 
E. Gabet, pers. comm.). In addition, 10 Be CRN ero- 
sion rate determinations were made throughout the 
catchment, both from bedrock exposures and 0 th - 
order (< 0.01 km 2 ) catchments, and from the mouth 
of the Khudi River where it joins the Marsyandi 
(Heimsath, in prep.). The measured CRN erosion 
rates were calculated using the same production rates 
and scaling factors as in the model, and have an av- 
erage error of 10%. Here we compare the statistical 
distribution of 56 10 Be CRN erosion rate determina- 
tions from bedrock exposures and 0 th -order catch- 
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Figure 8. Mean concentration of I0 Be across the simulation, showing the effect of initiation and cessation of landslides on 1 °Be concentration 
at several different rates of sediment detachment and landsliding. A) Mean I0 Be concentration through time for all model runs; initiation and 
cessation of landsliding indicated. B) Symbols highlighting variation in response time to changing boundary conditions, see text for discussion. 
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merits with predicted statistical distributions of CRN 
erosion rates determined from our simulation. Ad- 
ditionally, we calculate distributions of CRN erosion 
rates predicted at the mouth of the Khudi River, and 
compare these with the basin-wide CRN erosion rates 
determined for the Khudi catchment. 

The simulation of landsliding, cosmogenic produc- 
tion, and erosion was performed as described for the 
theoretical study in the San Antonio Creek watershed, 
with only a few modifications. First, the highest reso- 
lution DEM available for Nepal has roughly 90-m cell 
spacing, as opposed to the 30-m spacing available for 
San Antonio Creek. The model cell size was adjusted 
accordingly, and the smallest landslide populated in 
the model was also adjusted to have a radius of 90 
meters. Second, available erosion rate data from ther- 
mochronologic and sediment-load studies were used 
to establish a range of landsliding rates to input to the 
model of the Khudi catchment. A sediment detach- 
ment rate of 0.15 mm/yr was assumed for the initial 
model runs, and landsliding rates of 2.85, 3.35, and 
3.85 mm/yr (for total erosion rates of 3, 3.5, and 4 
mm/yr) were selected. The best-fit run was then re- 
analyzed with varying sediment detachment rates to 
derive a statistical distribution of bedrock and small- 
order catchment CRN erosion rates for comparison 
with existing CRN data. Finally, the size of the Oth- 
order catchments is of order the size of the model 
cell spacing, so these small catchments were treated 
as points within this model. 

3.4.1. Basin-wide CRN erosion rate in the Khndi 
catchment 

The distribution of basin-wide CRN-erosion rates 
derived from the model were compared with a 10 Be- 
derived erosion rate measured on sediment deposited 
at the mouth of the Khudi catchment. The measured 
CRN erosion rate of 3.4 mm/yr matches the median 
of the CRN erosion rates for the model run with total 
erosion rate of 3.5 mm/yr (0.15 mm/yr sediment de- 
tachment and 3.35 mm/yr landslides; Fig. 9A). This 
result supports the theoretical determination that at 
high landslide erosion rates, the effects of landslid- 
ing on CRN erosion rates are spatially averaged over 
large catchments, and that CRN methods are an effec- 
tive means of assessing average erosion rates. Further, 
the median values of each of these three model runs 
does not overlap the 1 st through 3 rd quartiles of any 


other run (Fig. 9A), indicating that the medians of 
these model runs statistically differ [23], The results 
of this simulation, in concert with the actual CRN ero- 
sion rate, allow a determination of the erosion rate in 
the Khudi catchment of 3 ± 0.5 mm/yr, tightening the 
constraints on the erosion rate as derived from ther- 
mochronologic and sediment-load data. 

3.4.2. Bedrock and small-order catchment CRN 
erosion rates in the Khudi catchment 

Using the basin-wide total erosion rate of 3.5 
mm/yr, four additional model runs were executed at 
increasing rates of sediment detachment (0.01, 0.05, 
0.10, and 0.15 mm/yr) to compare the simulated dis- 
tributions of CRN erosion rates with the distribution 
of CRN erosion rates determined from 56 10 Be ero- 
sion rates from bedrock samples and 0 th -order catch- 
ments (Heimsath, in prep.; Fig. 9B). In the model 
space, the same 56 locations that were actually sam- 
pled were used to determine a theoretical distribution 
of bedrock erosion rates. The results of this compar- 
ison indicate that, as hypothesized from the theoreti- 
cal results, sampling bedrock and small area drainage 
basins for CRN concentrations in regions with sig- 
nificant erosion by bedrock landsliding is not an ef- 
fective approach to determining basin-wide erosion 
rates. The data set analyzed in the Khudi drainage is 
one of the larger CRN erosion rate studies undertaken, 
yet the number of samples collected, and the percent- 
age of the landscape that they represent, is inadequate 
to provide a meaningful representation of basin-wide 
erosional processes. Although thermochronometric, 
cosmogenic, and sedimentologic techniques all indi- 
cate erosion rates in the Khudi of 2-5 mm/yr, the me- 
dian CRN erosion rate from the Khudi bedrock sam- 
ples is ~ 0.1 mm/yr. while the median of the simula- 
tion runs is ~ 0.25 mm/yr (Fig. 9B). It is also interest- 
ing to note that the median erosion rates from the four 
model runs are very similar, despite an order of mag- 
nitude difference in sediment detachment rate among 
them, and they differ significantly from the median 
erosion rate determined on the actual bedrock sam- 
ples. This difference can likely be ascribed to two 
sources, both due to the fact that the actual bedrock 
samples from the Khudi catchment almost all came 
from ridge crests, or small draws that head at ridge 
crests. The first source of the observed difference 
is that landslides in the model space are placed ran- 
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Figure 9. A) Distribution of model basin-wide CRN-derived erosion rates from the mouth of the Khudi drainage for erosion rates of 3.0, 3.5, 
and 4.0 mm/yr. The thick gray line across the center of the figure is the measured 10 Be derived erosion rate of 3.4 mm/yr from the mouth 
of the Khudi River (A. Heimsath, unpub. data). B) Comparison of simulated distributions of CRN-derived erosion rates from 56 sample 
locations (white boxes) where actual CRN-erosion rates were measured (gray box: A. Heimsath. unpub. data) for bedrock samples in O th -order 
catchments in the Khudi drainage. Median CRN erosion rates for both the simulated data and the actual CRN data are an order of magnitude 
less than the basin-wide erosion rate (A). 


domly, such that a ridge crest is as likely to be af- 
fected by landslides as a hillslope. In reality, this is 
likely to not be the case. Landslide mapping indicates 
that although landslides are preferentially located in 
areas of steep slopes near the headwaters of drainage 
systems, such slides rarely breach the drainage di- 
vide and lower the interfluve (e.g. [24]). Lowering 
of drainage divides most plausibly occurs during in- 
frequent events when the topographic slope below' the 
ridge crest has been over steepened by repeated land- 
sliding, and is thus this process is likely to be highly 
undersampled. The second effect is sampling bias. 
The CRN samples were specifically collected from 
bedrock outcrops along the crest that were judged to 
not have been recently affected by landsliding. As 
such, the spread of CRN ages would be expected to 
be smaller than a truly random set of ridge crest sam- 
ples, and the median erosion rate would, as a result, 
be lower. This difference notwithstanding, the results 
of this comparison still serve to emphasize the impor- 
tance of considering the effects of spatial and tem- 
poral averaging in collecting material for CRN deter- 
minations of erosion rates in actively landsliding re- 


gions. 

4. Discussion and Conclusions 

We have presented a numerical simulation for mod 
eling the production, decay, and distribution of cos- 
mogenic nuclides on a landscape, and their removal 
through sediment detachment and mass wasting pro- 
cesses. Although this simulation was developed to 
model the effects of landsliding on CRN-derived ero- 
sion rates, the framework of the simulation could 
be adapted to model any number of factors that ef- 
fect erosion rates derived from CRNs in sediments, 
including spatial variations in lithology and mineral 
content, ice cover, annual snowfall, and recent effects 
of glaciation. In the past, careful researchers trying 
to exploit CRNs to obtain erosion rates have com- 
monly restricted their sampling to small, unglaciated 
catchments with uniform lithologies and slow' ero- 
sion rates. Variations in CRN production due to to- 
pographic shielding, slope, lithology, snow and ice 
cover, and glacial history have generally been ig- 
nored. With the ability to numerically predict the 
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effects of such variations, studies can be expanded 
into larger catchments with higher erosion rates, di- 
verse sediment-production processes, and spatial het- 
erogeneities in CRN concentrations and events that 
can re-set the cosmogenic clock. 

Our landsliding model is underpinned by the ob- 
served frequency-magnitude relationships of land- 
slides in two mountain ranges. These happen to 
yield nearly identical exponential scaling factors (0 = 
1.1) that imply that large, infrequent landslides dom- 
inate the total sediment flux. The observation that 
implementation of our frequency-magnitude based 
landsliding module yields volumetric rates that os- 
cillate around the expected value suggests that the 
model succeeds in mimicking a natural process. This 
frequency-magnitude relationship and the appropriate 
scaling factor need to be verified in other mountain 
ranges before they are routinely applied. In addition, 
we make the implicit assumption that all of the sedi- 
ment generated by a landslide is delivered and homo- 
geneously mixed with sediment detachment-derived 
sediment within a single model time step (100 years). 
We make no attempt to model sediment storage on 
hillslopes or within fluvial systems, or to model parti- 
cle size fractionation between sediment detachment- 
derived and landslide-derived debris. 

Our specific evaluation of the effects of bedrock 
landsliding on erosion rates derived from CRNs us- 
ing this simulation yielded several results that have 
been previously described, particularly that stochas- 
tic processes, dominated by large rare events, are 
difficult to measure using basin-averaging sediment- 
sampling techniques, because the large, rare events, 
will often not be represented in the sample popula- 
tion (e.g. [?]). However, our simulations also indi- 
cate that the median CRN-determined erosion rates 
are representative of the volumetric erosion rates de- 
rived from the same catchments. This suggests that 
sampling multiple, similarly sized catchments, even 
in active, landslide-dominated mountain belts, offers 
a significant likelihood of yielding several samples 
with consistent CRN erosion rates. These rates are 
likely to be representative of the recent erosion within 
those catchments, although it must be recognized that 
such results will typically be lower than the long- 
term average that samples the large, rare events. As 
long as spatial variations in production can be ade- 
quately accounted for, as is done within our model- 


ing environment, then larger catchments will always 
yield a better approximation of long-term erosion rate 
in landslide-dominated terrains than smaller catch- 
ments. Although a specific relationship between the 
catchment size necessary to spatially average CRN 
samples and erosion rate is difficult to derive, based 
on our modeling results a general rule of thumb ap- 
pears to be 

Aw 9 = ^ (7) 

where A avg is the area needed to average the vari- 
ability in CRN concentration, and E is the estimated 
erosion rate over the catchment. 

Sampling from bedrock outcrops to measure aver- 
age erosion rates in landslide-dominated catchments, 
however, is unlikely to be a useful exercise. At the low 
end of landslide erosion rates, such samples will faith- 
fully yield the sediment detachment rate on the land- 
scape, but at increased rates of landslide erosion, such 
samples will only yield a very rough upper bound on 
the sediment detachment rate, and, more likely, will 
be uninterpretable in the context of the spatially aver- 
aged erosion rate over the study area. 

Finally, the response time of 10 Be concentrations 
over the landscape to changes in rates of erosional 
processes is thousands to tens of thousands of years. 
Regions that have undergone recent changes in rates 
of erosion may yield CRN-derived erosion rates that 
reflect some intermediate rate between the previous 
and current erosion rates during the re-establishment 
of ln Be equilibrium over the landscape. 
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Appendix A 
Governing Equations 

This Appendix details the derivations of three spe- 
cific portions of the landslide model, calculating the 
rate of landsliding, k, given a mean landslide erosion 


16 


Niemi, Oskin, Burbank and Heimsath 


rate, Ei s \ population of landslides in model space, and 
determining mean cosmogenic concentration per vol- 
ume of material removed from the model. 

Rate of Landsliding 

We assume that landslides follow a power law fre- 
quency magnitude distribution [7,14]. For such a dis- 
tribution, the cumulative frequency of landslides can 
be written as 

tia>a. = * (A a /A T )~ e A r , (A8) 

where tia>a. is the number of landslides greater than 
area A s that occur in a given year, A r is a specified 
reference area, k is the rate of landsliding, and 8 is 
the power-law exponent of the frequency-magnitude 
distribution [7, see eqn. 1], If A r is taken to be 1 km 2 , 
then from Eqn. A8, the total number of landslides per 
year, n, over the reference area can be written as 

n T = f nA~ p dA, (A9) 

“ Amin 

= K | A-efr", (A 10) 

= *{KL-Ki t)- (AH) 

where A min is the area of the smallest landslide to oc- 
cur in A r (in the case of the simulation, A min is equal 
to the model cell size) and A max is the largest slide 
area, constrained by local topographic relief. Since 
A m i n and A max can be constrained either empiri- 
cally or based on physical characteristics of the modei 
space, and /? has been shown to be ~ 1 over a large 
range of landslide erosion rates [7.14,15], the volume 
of material removed by landsliding must be controlled 
by the rate of landsliding, k. Because we would prefer 
to prescribe a rate of erosion due to landslides, Ei s for 
the model, we need to solve for k in terms of Ei s . To 
do this, we start by determining the number of slides, 
tia „, of a given area, A s per year over the reference 
area, A r : 

tia , = (A 12) 

Given this, the volume of erosion in any year due to 
landslides of area A s is equal to the number of land- 
slides of area A s multiplied by the volume of land- 
slides of area A„, 

V A , = tia. • V s (A 13) 

= nPA-'-e-V.. (A 14) 


In this simulation, we have opted to model the land- 
slides with parabolic cross sections and a linear re- 
lationship between maximum slide depth and width 
[25,7]. The scaling between landslide area and depth 
is defined by a scaling factor, £, with an empirically 
determined value of ~0.05 [7]. Thus, the depth of 
the landslide, d, at any radial distance from the land- 
slide center is a function of the maximum slide depth, 
the radial distance from the landslide center, r, and a 
constant, C, 

d = d max — Cr 2 . (A15) 

We can solve for C at the outer edge of the landslide, 
r m «i, where d = 0, 


For a slide of a given area A s , then, the maximum 
landslide depth, d max and radius, r max , are given by 


dmax — £ x/ As 


Substitute these values into Eqn. A16 to solve for C 
in terms of area, A s : 


and substitute Eqns. A17 and A19 into Eqn. A15 to 
solve for d as a function of A s , 

d = £-/a 7 7==r 2 . (A20) 

VAs 

The volume, V s , of a landslide of area A s can then be 
calculated by integrating over cylindrical shells from 
V — 0 tO T ma x , 


r / SyJ~A s - 
J 0 


=r 2 • rdr (A21) 


0 | 2 ^4 

= 2 tt £ r 7=r 

I 2 4VC47 0 

9 (aT ^ /2 \ 

= 2n£ {-2^--^r) 

= Ui'* 

2 8 
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Eqn. A24 can be substituted back into Eqn. A 14 to 
yield the volume of material, Va., removed from the 
reference area, A T per year by slides of area A s \ 

V A . = \nep AW*-®. (A25) 

The total volume of erosion per year from reference 
area A r , then, is the sum of the erosion due to land- 
slides of all sizes that occur in A r , from A m i n to 
Amax > 




1 r^max 

Ike/ 3 [ A\] /2 - B) dA (A26) 

2 JA mi „ 

«A27, 


The total erosion rate per year due to landslides, then, 
can be determined by dividing the volume of material 
eroded per year by landslides over A r by the area of 
A r , 



(A28) 


and, since we have defined A r to be a unit area (1 
km 2 ). 


Eu = V r (A29) 

Thus, Eqn. A27 can be solved for k in terms of £) s to 
yield 


* = -/ 'I, , )V (A30) 

Substituting Eqn. A30 into Eqn. All yields the to- 
tal number of landslides, n r , per year over the refer- 
ence area, A r in terms of known or prescribed values 

Amin, Am ax , &■> Und El s ', 


Tl r 


E u (3-2/3) (A-^-A-f n ) 


(A31) 


The total number of landslides, n; s , then, in the simu- 
lation for a given time step, t, over the entire simula- 
tion is 


n ls = t x n T x 4=2. (A32) 

A r 

where A S j m is the area, in km 2 , of the simulation. 


Landslide Population and Distribution 

With the number of landslides per time step deter- 
mined, the simulation is populated. The position each 
landslide is specified by a randomly generated x,y co- 
ordinate pair. The size of a given landslide is derived 
from the landslide frequency-magnitude relationship 
(Eqn. A8). The probability of a landslide with area 
A„ occurring is 

Pa. = « A~ 0 (A33) 

therefore, randomly generated numbers mapped lin- 
early onto the range Pa„ 4 „ to P Ama , can be used 
to create a population of nu landslides that fit the 
frequency-magnitude distribution nA~ 0 . 

Cosmogenic Nuclide Concentration in Eroded Ma- 
terial 

For each model step, the total depth of material re- 
moved from a given cell, D, is the sum of material 
eroded by sediment detachment, E s and the material 
eroded by landslides, Ei s . Using the surface con- 
centration of a cosmogenic nuclide at the end of the 
model step, A r ^°Be, the average nuclide concentration 
in the volume of eroded material, AT^Be can be cal- 
culated by integrating over the total depth of eroded 
material the concentration of the nuclide as a function 
of e-folding depth, A, 


AT)°Be • p , 

[ e~ z ^ A) dz 

(A34) 

O 

-D 


Nj°Be ■ p ■ 

A L e -z(p/A)|° 

(A35) 


p \ 1 -D 


NfBe • A ( 

r e D(p/A) _ 

(A36) 
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Abstract 


The Kyrgyz Range, the northernmost portion of the Kyrgyzstan Tien Shan, displays topographic 
evidence for lateral propagation of surface uplift and exhumation. The highest and most deeply 
dissected segment lies in the center of the range. To the east, topography and relief decrease, and 
preserved remnants of a Cretaceous regional erosion surface imply minimal amounts of bedrock 
exhumation. The timing of exhumation of range segments defines the lateral propogation rate of the 
range-bounding reverse fault and quantifies the time and erosion depth needed to transform a 
mountain range from a juvenile to a mature morphology. New apatite fission-track (AFT) data from 
three transects from the eastern Kyrgyz Range, combined with published AFT data, demonstrate 
that the range has propagated over 110 km eastwards over the last 7-1 1 Myr. Based on the thermal 
and topographic evolutionary history, we present a model for a time-varying exhumation rate 
driven by rock uplift and changes in erodability and the time scale of geomorphic adjustment to 
surface uplift. Easily eroded, Cenozoic sedimentary rocks overlying resistant basement control 
early, rapid exhumation and slow surface uplift rates. As increasing amounts of resistant basement 
are exposed, exhumation rates decrease while surface uplift rates are sustained or increase, thereby 
growing topography. As the range becomes high enough to cause ice accumulation and develop 
steep river valleys, fluvial and glacial erosion become more powerful and exhumation rates once 
again increase. Independently determined range-normal shortening rates have also varied over time, 
suggesting a feedback between erosional efficiency and shortening rate. 

Introduction 

Growth of a contractional mountain range is driven by an evolving relationship between rock 
uplift, surface uplift, and exhumation. Surface uplift is defined as rock uplift minus exhumation 
[. England and Molnar, 1990]. In the early phase of orogenesis, rock uplift must outpace 
exhumation; in a steady-state orogen, these reach an equilibrium, while in the final, destructive 
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phase of an orogen, exhumation dominates [e.g., Willett and Brandon, 2002], Temporal changes in 
key controlling factors, such as fault geometry, shortening rate, and erosion rate will influence the 
surface uplift history [e.g., Burbank et al., 1996; Abbott et al., 1997; Stock atui Montgomery, 1999; 
Whipple and Tucker, 1999]. Erosion rate may be controlled by factors such as channel gradients, 
surface slopes, relief, precipitation, glacial extent, and rock resistance [e.g., Ahnert, 1970; Howard, 
1994; Brozovic et al.. 1997; Hallet et al. 1996; Schmidt and Montgomery , 1996; Sklar and Dietrich, 
2001; Whipple, 2004]. Perhaps because complex interactions among these factors undoubtedly 
occur, the detailed evolution and interdependence of these factors is rarely delineated. 

Reverse-fault-bounded mountain ranges propagating into a foreland basin commonly initiate 
either as a single, localized structure which gradually lengthens along strike with increasing amount 
of shortening, or as several fault segments which eventually coalesce or overlap [ Dawers , 1993; 
Cartwright et al., 1995], In either case, the active range-front likely lengthens as displacement 
accumulates on the range-bounding faults. After several million years of displacement, the 
original, relatively small structures may be impossible to discern. However, if a reference horizon 
along the trend of the range exists, range growth and exhumation can be placed into a topographic 
reference frame by combining low-temperature thermochronologic data with structural geology and 
geomorphic analysis. When ranges grow through lateral and vertical propagation, a space-for-time 
substitution can illuminate the long-term spatial and temporal distribution of exhumation and 
surface uplift and can permit reconstruction of progressive changes in the balance between rock 
uplift and exhumation [ Burbank et al., 1999]. Space-for-time substitutions are most reliable when 
time constraints exist for the interval of range propagation. Such along-strike temporal control, 
however, is commonly lacking in most studies of fold-and-thrust belts. Only when reliable ages can 
be assigned to various stages of range growth can the validity of the substitution be assessed. We 
report here a time-calibrated example of range propagation from the Kyrgyz Tien Shan. 
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The Kyrgyz range of the northern Tien Shan (Figure 1) provides an example of a reverse-fault 
bounded mountain range with topographic evidence for progressive lateral propagation of surface 
uplift and exhumation. The highest and most deeply dissected sector of the range lies in its central 
portion, south of Bishkek, in the region of the Ala Archa River. Glaciers mantling peaks up to 4800 
m have deeply dissected the granitic range in this region, and reset apatite fission-track and U- 
Th/He thermochronometers indicate >5 km of exhumation at river level [ Bullen et al., 2003]. To 
the east and west, both topography and relief decrease. Moreover, preservation of remnants of a 
Cretaceous regional erosion surface [ Trofimov et al., 1976] implies minimal amounts of bedrock 
exhumation. If the timing of exhumation of segments of the range outward from the center can be 
defined, the lateral propagation rate of the range can be estimated and the depth of erosion and time 
needed to transform a mountain range from a juvenile to a mature morphology can be 
reconstructed. 

Here we present new apatite fission-track (AFT) data from three transects from the eastern half 
of the Kyrgyz Range that, when combined with published AFT data, demonstrate that the range has 
propagated over 110 km from the presently highest region towards the east over the last 7-1 1 Myr. 
When synthesized with structural data and analysis of recently produced digital topography, we can 
reconstruct the thermal and topographic evolutionary history of the Kyrgyz range. We present a 
model for a time-varying exhumation rate driven by the interplay of rock erodability, surface 
processes, and shortening rate. We propose that easily eroded Cenozoic sedimentary rocks that 
overlie resistant pre-Tertiary bedrock control early, rapid exhumation and slow surface uplift rates. 
As increasing amounts of resistant bedrock are exposed, exhumation and shortening rates decrease 
while surface uplift can persist, thereby growing topography. This decrease in shortening rate may 
be counterbalanced by basinward-propagation of the deformation front. As the range becomes high 
enough to cause ice accumulation and orographically enhanced precipitation, fluvial and glacial 
erosion become more effective and exhumation rates once again increase. A balance between 
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erosion and rock uplift occurs through establishment of mature drainage networks into the range 
that connect high-elevation glaciated areas to deeply incised rivers. Based on analysis of the digital 
topography, we argue that the central part of the range has achieved this balance and represents a 
topographic steady state. 

Geologic history 

The Tien Shan records a complex Paleozoic histoiy of island arc accretion [Burtman 1975; 
Carroll et al., 2001; Bazhenov et al., 2003] followed by Permian strike-slip deformation [Burtman 
1975; Bazhenov et al., 1999]. Episodes of intracontinental deformation driven by distal plate 
margin tectonism occurred during the Early-Middle Jurassic and the Late Jurassic-Cretaceous, 
documented by foreland-basin formation to the north and south of the Tien Shan, as well as in a 
prominent transtensional basin cross-cutting the range [Hendrix et al., 1992; Sohel, 1999], Apatite 
fission-track data suggest pulses of exhumation during the Permian and the Jurassic [Sobel and 
Dumitru, 1997; Bullen et al., 2001; Dumitru et al., 2001]; these are likely correlated with episodes 
of deformation within the Tien Shan and deposition within the Tarim basin. 

A widespread regional erosion surface formed within the central Kyrgyz Tien Shan in the late 
Mesozoic [Trofimov et al., 1976; Makarov, 1977; Chediya, 1986; see Burbank et al., 1999 for 
photographs]. This surface is unconformably overlain by the Paleocene-Eocene Suluterek 
formation (also called the Kokturpak formation), containing calcareous sandstone, dolomite, 
gypsum, and an arid spore-pollen assemblage [Chediya et al., 1973; Fortuna et al., 1994] as well as 
localized Eocene basalt flows [Krylov, I960]. The formation is typically ca 150 m thick and is 
truncated by an erosional unconformity separating it from upper Tertiary sediments [Chediya et al., 
1973; Fortuna et al., 1994]; drilling in the Chu basin reveals that the thickness locally reaches 635 
m [Burg et al., 2004], Overall, this Mesozoic to Tertiary interval of erosion and minor deposition 
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spans ca. 100 Myr period of tectonic quiescence, during which a stable thermal regime in the upper 
crust was established [Bullcn et al., 2001], 

The Cenozoic exhumation history of the Tien Shan has been used to study the evolution of the 
range. Fission-track cooling ages and deposition of apparently syntectonic conglomerates in the 
adjacent Tarim basin suggest that shortening commenced around the Oligocene - Miocene 
boundary [Hendrix et al., 1994; Sobel and Dumitru , 1997; Yin et al., 1998], This intracontinental 
deformation is driven by the Eocene - present collision of India with Asia. Deformation within the 
Tien Shan between 73-80°E longitude (Fig. 1) appears to have begun along the southern side of the 
range adjacent to the Tarim Basin at ca. 26 Ma and then propagated northwards across the 
individual ranges of the Tien Shan [Sobel et al., 2000; Dumitru et al., 2001]. Deformation has not 
propagated monotonically northward; at least the last several million years of this history have been 
marked by deformation throughout the entire range [Thompson et al., 2002]. At present, seismicity 
is distributed throughout the width of the orogen [Bune and Gorshkov, 1980] and geodetic studies 
document a continuous gradient of shortening between Kashgar and Bishkek [Abdrakhmatov et al., 
1996; Reigber et al., 2001], Exhumation of the Kyrgyz range on the northern margin of the Tien 
Shan began in the region of Ala Archa (Fig. 2) at ~1 1 Ma [Bullen et al., 2001], Deformation is 
dominantly north-vergent thrusting, with a minor sinistral transpressional component [e.g., 
Cobbold etal., 1996; Mikolaichuk, 2000; Thompson et al., 2002] (Fig. 2). 

The Oligocene-Miocene Shamsi Formation, exposed in the southern margin of the Chu Basin, 
was deposited in a foreland basin that pre-dated growth of the Kyrgyz Range [Bullen et al., 2001]. 
The thickness of the pre-1 1 Ma sediment in the Chu basin is ca. 1 km. Although it is difficult to 
constrain the amount of Cenozoic sediment that formerly overlay the Kyrgyz range, it likely 
exceeded the amount in the Chu basin because of flexural subsidence of the Kyrgyz Range in 
advance of the northward propagating Tien Shan uplift. This pre-1 1 Ma burial contributed to the 
maximum temperature experienced by basement samples during the Neogene. During subsequent 
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exhumation, the weakly lithified Neogene sediment could be more readily eroded than the 
underlying Paleozoic strata. The Chu basin is deepest in the region of Bishkek and becomes 
shallower to the east, west, and north [ Abdrakhmatov et al., 2001] (Figure 2), suggesting that the 
magnitude of flexural subsidence and hence the amount of thrusting in the Kyrgyz range is greatest 
in the central section of the range. Bullen et al. [2001; 2003] present a detailed basin and tectonic 
analysis based on the stratigraphy and magnetostratigraphy of the last ca. 9 Ma of the Chu basin fill 
combined with structural and apatite U-Th/He and fission-track thermochronologic studies of the 
central Kyrgyz range. This work demonstrates that the range was rapidly exhumed between 1 1 and 
10 Ma along a north-vergent thrust; deformation and cooling rates decreased significantly for the 
next 7 Myr until 3 Ma, when rates increased again. The deformation front propagated slightly 
between 8 and 3 Ma, coinciding with coarsening-upwards deposition in the basin. During the last 3 
Myr, sediment accumulation rates reached a maximum, while the Kyrgyz range was being rapidly 
exhumed; deformation propagated farther into the foreland, forming a piggyback basin. 

Geomorphology and Neotectonics 

GPS measurements indicate ~20 mm/yr of modem shortening across the central Tien Shan: 
nearly half of the total convergence rate between India and Asia [Abdrakmatov et al., 1996; 
Reigber et al., 2001], This shortening is distributed across a 400-km-wide belt of subparallel 
ranges and intramontane basins that are separated by active reverse faults [Thompson et al., 2002]. 
The Kyrgyz Range forms the northernmost topographic crest within the central Tien Shan. Young 
thrust-fault scarps and seismicity along the northern margin of the range [Chediya et al., 1998; 
Thompson et al., 2002] attest to ongoing shortening. These ranges form an orographic barrier, with 
precipitation focused on the northern flanks [Aizen et al., 1995]. 

From a structural perspective, a propagating range would be expected to exhibit increasing 
amounts of rock uplift and structural relief as a function of distance from the tip toward the center 
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of the range. From a geomorphological perspective, the temporal lag between surface uplift and 
erosional response [Kooi and Beaumont, 1996; Willett, 1999] suggests that relief, peak height, 
hillslope angles, and mean elevation should all initially increase as the range grows. As a range 
approaches or attains a topographic steady state, relief, mean elevation, and hillslope angles should 
stabilize. Hillslope angles, topographic relief, and hypsometry were derived from the 3-arcsec 
DEM (SRTM, 2003) resampled into Universal Transverse Mercator projection with 70 m pixels. 
Hillslope angle, mean and maximum elevation of the Kyrgyz Range all attain average values within 
50 km of its eastern tip, supporting the contention that range-scale topographic steady state may 
prevail over the majority of the range (Fig. lc). 

A view of the early stages of geomorphic evolution of the Kyrgyz range is provided by 
observations of uplifted and deformed remnants of the pre-Cenozoic erosion surface on its 
southeastern slope (Fig. 3). This exhumed unconformity surface is recognized in the field and on 
remote sensing images as concordant, uniformly tilted regions etched by a distinctive dendritic 
network of shallow bedrock channels [Os kin and Burbank, Alpine landscape evolution dominated 
by cirque retreat, submitted manuscript to Geology, hereafter referred to as Oskin and Burbank, 
submitted manuscript ]. On the southern flank of the Kyrgyz Range, southward tilting of the erosion 
surface remnants increases from east to west, consistent with a model of increasing displacement 
and range-scale limb rotation above a listric thrust fault [Erslev, 1986], such as is proposed to 
underlie the Kyrgyz Range [Abdrakhmatov et al., 2001]. Field observations at the southern foot of 
the range and at other nearby outcrops of the erosion surface indicate that it is exhumed from 
beneath easily eroded Cenozoic continental sedimentary rocks. Remnants of the erosion surface can 
be traced from the foot of the range to the range crest in the easternmost south-facing slope of the 
Kyrgyz Range. Outcrops of the erosion surface diminish westward in response to progressively 
greater range uplift above the Pleistocene glacial equilibrium line altitude (ELA), and 
corresponding incision of glaciated valleys into bedrock. North-facing glaciated valleys also 
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expand laterally via cirque retreat into low-relief south-facing erosion surface remnants [Oskin and 
Burbank, submitted manuscript]. This process expands the area covered by north-flowing glaciers 
(Fig. 3) and moves the drainage divide relatively southward. 

Morphometric analyses of drainage basins that descend from the crest of the Kyrgyz Range link 
catchment-scale erosion and relief generation to range-scale topographic development. Most of our 
morphometric analyses focused on the north-facing catchments that reach the range crest. Basins 
were each derived with standard hydrologic routing functions and hypsometry and hillslope angles 
(Figure 4) were measured directly from the topography and its derivative, respectively. Relief 
measurements, calculated as the difference between two elevations, are inherently subject to bias 
from the horizontal length scale between elevation measurements. Rather than measure mean relief 
over a fixed distance, which would primarily reflect mean hillslope angle, an alternative relief 
measure was devised to convey the degree of overall range dissection. This relief measurement, 
termed here as the internal relief, is defined as the maximum elevation difference among all 
elevation points equidistant from the basin outlet, e.g., the lowest and highest elevation among all 
the points that are 5 km from the outlet. Typical values of internal relief climb monotonically 
upstream from the basin outlet, hover around a maximum level in the middle basin reach, and 
descend towards the divide (Fig. 5). In practice, internal relief correlates to the product of basin 
area and mean hillslope angle, and therefore provides a measure of the ridge-to-valley relief that 
reveals how extensively glaciers and rivers have incised below adjacent peaks. Because extreme 
values of relief and elevation within Kyrgyz Range catchments tend to bias values toward isolated 
high peaks that occupy only a small fraction of the landscape, the upper quartile of internal relief 
and the lower, middle, and upper quartiles of elevation (as defined by the hypsometry) provide a 
better characterization of key topographic attributes (Fig. 4). 

The distributions of hillslope angle, hypsometry, and internal relief from east to west within the 
Kyrgyz Range reveal a pattern of initial surface uplift followed by dissection and then stabilization 
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of range-scale morphology. At the easternmost end of the range (surface uplift zone on Fig. 4), the 
upper quartiie of internal relief is < 600 m, the mean elevation is <3 1 00 m, peak heights are <4000 
m, and the range of hypsometry (measured by difference between the 1 st and 3 rd quartiles of the 
hypsometry) is <500 m. Low relief and a limited range of altitude indicate only limited dissection 
of this high mountain range and support predominance of surface uplift over erosion here. Moving 
westward, peak elevations, internal relief, and the range of hypsometry all continue to increase 
toward the center of the Kyrgyz Range (transition zone on Fig. 4), whereas the mean elevation 
remains nearly constant across this 80-km span. The increase in the elevation of the highest peaks 
suggests that total rock uplift increases toward the west Similarly, the growth of internal relief and 
hypsometric range indicate increasing dissection across this tract. In this zone, adjacent basins 
compete to attain necessary size and gradients in order to balance rock uplift. Overall stability of 
morphologic indices along the central part of the range (steady morphology zone, Fig. 4) suggests 
that this region may have reached an equilibrium topographic form in which erosion by balances 
tectonic uplift. 

Fission-track methodology 

Apatite fission-track data from exhumed basement rocks often yield distinctive age-elevation 
patterns that can be used to infer the low-temperature exhumation history of a range. Samples may 
have been exhumed from sufficient depth during the most recent exhumation event that the 
maximum temperature, T max , exceeded the total annealing temperature (Figure 6a). In this case, the 
fission-track clock was reset to zero, and fission-track data record information on the time- 
temperature cooling path of the sample as it cooled through the partial annealing zone (PAZ) 
during exhumation [e.g. Green et al., 1989a, 1989b]. For apatites that are cooled moderately 
rapidly (10°C/Myr) and which have kinetic parameters similar to Durango apatite, the total 
annealing temperature is ~120°C [Donelick et al., 1999; Ketcham et al., 1999]. Samples in the rock 
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column above the total annealing isotherm experienced lower T max prior to exhumation; these 
samples resided in the (now exhumed) PAZ for some period of time and therefore have ages 
reflecting the penultimate cooling event, (strongly) modified by partial annealing [e.g., Fitzgerald 
et al., 1995] (Figure 6a). The age of the transition between these two suites of samples is 
interpreted to represent the onset of rapid exhumation. 

The apatites analyzed in this study often have variable kinetic properties, as documented by 
both etch pit diameter (Dpar) and Cl content. These two methods have been shown to yield equally 
useful data for assessing the kinetic properties of apatite [Donelick et al., 1999; Ketcliam et al., 
1999]. In particular, apatites distinguished by large etch-pit diameters and high Cl content (herein 
termed “more resistant”) are typified by higher annealing temperatures in contrast to less resistant 
apatites [e.g., Green et al., 1989b; Ketcham et al., 1999], Given these contrasting annealing 
temperatures a superficial interpretation of an age-elevation plot can yield false conclusions about 
both the onset and rate of exhumation. Although the total annealing temperature of different 
apatite types is variable, the low temperature annealing behavior of apatites appears to be similar 
[Ketcham et al., 1999], Therefore, the slope of the exhumed partial annealing zone on an age- 
elevation plot for different apatite types cannot be assumed to be parallel (Figure 6). To avoid this 
pitfall, we have plotted age-elevation curves for apatites with similar kinetic properties. The long 
period of tectonic quiescence prior to the onset of Neogene exhumation implies a stable thermal 
structure in the upper crust [Bullen et al., 2001], Therefore, the slope of the age-elevation curve 
within the exhumed PAZ for kinetically similar apatites which have experienced similar thermal 
histories should be about the same for all of our profiles. This assumes that only a small, similar 
amount of horizontal-axis rotation has been experienced at each of the profiles located in the 
hanging wall of a major thrust. A contrasting age-elevation slope will pertain for apatites with 
differing kinetic characteristics (see Figure 6). 
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Sample preparation and analytical details are presented in Table 1. The young apatite samples 
yielded very few horizontal confined track-length measurements; up to 100 track-lengths were 
measured from older samples. For age determinations, 15 to 31 grains per sample were selected at 
random and dated; 1 sample yielded only 4 countable grains. Following convention, all statistical 
uncertainties on pooled ages and mean track lengths are quoted at the ±lo level, but ±2o 
uncertainties are taken into account for geologic interpretation. To assess the kinetic properties of 
apatite, four Dpar measurements were averaged from each dated crystal and from each crystal 
which yielded a confined track length, provided that sufficient etch pits were present. Dpar values 
are operator and etchant-dependent [Sobel et. al., 2004], Therefore, seven samples were also 
analyzed with a CAMECA SX-50 or SX-100 electron microprobe in order to determine Cl content 
(Table 1). Every crystal with a single grain age or a confined track-length measurement was 
probed. The microprobe data were used to calibrate the Dpar measurements. 

In the context of age-elevation data, we utilize the AFT data from 3 transects to reconstruct the 
onset and rate of exhumation for each transect. Subsequently, Cenozoic burial and exhumation 
histories are evaluated using thermal models of higher elevation samples. 

Results 

Shamsi 

Samples were collected along a tributary of the Shamsi river valley between 2250 and 3770 m 
along a transect that traverses a continuous sequence of Carboniferous quartzite from the base of 
the valley up to the crest of the east flanking ridge (Fig. 3). 

The six samples yielded Miocene to Jurassic ages, generally increasing with elevation (Table 1; 
Figure 7. Four of them define a readily interpretable succession that captures the base of an 
exhumed partial annealing zone. These samples also illustrate the complications introduced by 
apparent mixtures of apatites with different annealing temperatures within the same bedrock 
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sample. With an age of 6.3 ± 0.8 Ma, the basal, youngest sample (TS164) passes the % 2 test; all 
analyzed grains have similar kinetic properties as shown by Dpar (1.74 pm, SD 0.11) and 
microprobe (0.07 wt% Cl, SD 0.05) (Table 1; Figure DR1). The next sample (TS165) has a central 
age of 14 ± 7 Ma and fails the %1est. Eleven of the twelve countable grains form a young 
population, pass the x 2 test with a pooled age of 6.7 ± 1.8 Ma and a mean Dpar of 1.74 pm. A 
single grain has an age of 124 ± 44 Ma and a Dpar of 2.04 pm. However, microprobe data does not 
differentiate this grain; the average chlorine value of the sample is 0.12 wt% and the old grain has a 
value of 0.1 1 wt%. The third sample (TS166) passes the % 2 test and has a pooled age of 17 ± 2 Ma 
and a Dpar value of 1.74 pm. The highest elevation sample (TS170) fails the y 2 test and has a 
central age of 91 ±6 Ma, a Dpar value of 1.65 pm, and 0.12 wt%cl. However, one crystal has an 
anomalously high wt%cl of 0.45%; excluding this grain, the sample passes the % 2 test and has a 
pooled age of 88 ± 5 Ma and 0.07 wt%cl. For these latter two samples, there is no relationship 
between age and Dpar. Excluding the single old, high Dpar grain in sample TS165, these 4 
samples can be plotted together on an age-elevation plot (Figure 7). 

The lower two samples show rapid cooling, while the upper two show slow cooling within an 
exhumed PAZ. The inflection point between the two segments, at -2800 m and 7 Ma, defines the 
onset of rapid cooling. Assuming that the sample experienced some amount of heating due to 
Miocene sedimentation prior to rapid exhumation, the total annealing temperature corresponding to 
the inflection point of this low Cl path is 100-105°C. Based on the age and elevation differences of 
these two partially reset samples, the apparent exhumation rate for the exhumed PAZ is about 12 ± 
1 m/Myr (0.012 ± 0.001 km/Myr). This value will be used to define the apparent exhumation rate 
for the exhumed PAZ at Issyk Ata, because this section experienced a similar burial history (c.f, 
Fig. 6). As will be shown below, the corresponding rate at Boom Gorge is ca. 0.01 km/Myr, 
supporting this assumption. After 7 Ma. the rate at Shamsi accelerated to -1 km/Myr. 
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Two additional samples lie within the exhumed PAZ. Both pass the chi-squared test. TS167 
has an age of 67 ± 3 Ma, Dpar of 1.96 pm and 0.10 wt%Cl (0.08 wt%Cl excluding 5 grains with 
high values); TS169 has an age of 151 ±6 Ma, Dpar of 2.04 and 0.26 wt%Cl (Table 1; Figure 
DR1). Because TS169 is significantly older than the two bounding samples TS167 and TS170, 
plotting all of the samples within the exhumed PAZ together cannot yield a readily interpretable 
age-elevation plot. However, placing kinetically different samples on subparallel trends reveals 
that more resistant apatites have consistently older ages, representing older portions of a single 
exhumed PAZ. 

Issyk Ata 

Samples were collected from a transect along the Issyk Ata river valley between 1 840 and 3290 
m. (Fig. 4). Nearby, glaciated peaks reach elevations of 4500 m. The transect samples primarily 
Late Ordovician granite; the second sample above the base is an apatite-poor Riphean metavolcanic 
rock. 

Five of the six analyzed samples pass the chi-squared test, yielding ages ranging from 3.9 ± 0.7 
Ma to 6.9 ± 0.6 Ma (Table 1; Figure 7). Four of these samples, Mav38, TS158, TS159, and TS162 
appear to be monocompositional based on Dpar measurements and microprobe analysis. Dpar 
values are 1.65, 1.36, 1.87, and 1 .92,respectively; the former sample yields a wt%Cl of 0.02 (Table 
1; Figure DR1). Sample TS163 yielded only 4 countable grains; the resulting age is consistent with 
the other samples but of too low precision to warrant further attention. Sample TS 1 61 contains two 
components that each pass the chi-squared test. A young population of 25 grains yields an age of 
7.6 ± 1.9 Ma; an older component with 6 grains provides an age of 102 ± 11 Ma. The younger 
population has a Dpar of 1.72 pm and a wt%Cl of 0.02. Kinetic characteristics do not explain all of 
the grains in the older population: 5 of the 6 older grains yield higher values of 0.17 wt%Cl, 
whereas 1 grain has a value of 0.01 wt%Cl. Three of the old grains have large Dpar values of 2.43; 
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the other 3 fall within the cluster of young grains with an average value 1.83. When combined, the 
6 grains have an average value of 2.13 pm. 

The young, low resistant grains from these 6 samples define a steep linear trend on an age- 
elevation plot, indicating rapid cooling from below the base of an exhumed PAZ. The trend of the 
corresponding exhumed PAZ is assumed to have the same slope as the corresponding curve at 
Shamsi. Due to the absence of low resistance, partially reset samples, the position of this latter 
curve is poorly defined; hence, the onset of this rapid cooling can only be constrained as older than 
about 8 Ma. Using Fig. 6D as an analogy, the slope of the exhumed PAZ through the more resistant 
component of sample TS161 should have the same slope as the less resistant component; the two 
curves should be separated by a vertical distance corresponding to the difference in T a for the two 
components. Although the more resistant slope is poorly constrained, this analysis suggests that 
the onset of rapid cooling should be represented by an elevation close to TS161 and hence only 
slightly older than 8 Ma. 

Boom gorge 

Two samples were collected from the Boom gorge along the Chu River at 1375 and 1530 m; 
these lie ca. 1500 m below the local peaks (Fig. 4). Both samples are from topographically low 
positions in the footwall of the main range-bounding thrust, in contrast to the hanging wall sections 
sampled at Shamsi and Issyk Ata. Cenozoic sediments lying above the regional erosion surface are 
preserved in the center of the range, overthrust by Paleozoic strata. The same sedimentary sequence 
lies north and south of the range, overthrust from the south and north, respectively. The Cenozoic 
strata in this region do not exceed 1.5 km in thickness [ Trofimov et al., 1976], Sample TS84 was 
collected from a Permian granite just below the erosion surface; both the granite and the Cenozoic 
strata have been overthrust by Paleozoic units. Sample TS27 was collected from Devonian - Upper 
Carboniferous sandstone. 
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The two samples both pass the chi-squared test, with pooled ages of 128 ± 10 and 150 ± 8 Ma 
(Table 1; Figure 7). The older sample, TS84, has a Dpar of 1.70 pm (Table 1; Figure DR1). The 
younger sample, TS27, was prepared with slightly different etching conditions; therefore, Dpar was 
not measured. The samples yielded similar track-length data of 12.50 and 12.82 pm, respectively. 
The two samples lie on a cooling trend representing an exhumed PAZ; assuming that the two 
samples have the same kinetic characteristics, they experienced an apparent exhumation rate of 
0.01 km/Myr. Because sample TS27 was collected from a structurally deep position within the 
Boom gorge, the amount of exhumation in the gorge is clearly limited. 

Thermal modeling 

Track-length modeling cannot determine a unique thermal history; rather, it yields a range of 
solutions which are consistent with the observed data. Thermal modeling combining track length, 
single crystal ages and weight percent chlorine were performed using the AFTSolve program 
[Ketcham et al., 2000] and the annealing model of Ketcham et al. [1999], This program determines 
the best-fit models as well as good and acceptable fits. Modeling of partially annealed samples 
provides good constraints on the total annealing temperature and the maximum Miocene burial 
temperature; the minimum temperature prior to Late Cenozoic burial is less well constrained. With 
information on the thermal state of the crust, this information delineates the depth of the sample 
below the regional erosion surface and the Late Cenozoic sediments that formerly overlay the area. 
The constraints on the time of final exhumation of partially annealed samples from modeling are 
often less precise; however, this information may be available from lower elevation samples in the 
same transect. 

Converting thermal histories into exhumation paths requires information on the thermal state of 
the crust and, ideally, thermal conductivity. Borehole data suggest that the geothermal gradient 
within the study area is presently 25-30°C/km [Gubin, 1986]; herein, we use a value of 26°C 
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[Shvartzman, 1992]. Conductivity data are sparse. Exhumation rates are more difficult to calculate, 
given that one must also consider advection of isotherms due to thrusting [e.g.. Brown and 
Summerfield, 1997; Mancktelow and Grasemann, 1997] and perturbations of isotherms due to 
topographic effects [e.g., Stiiwe et al., 1994], In this study, neither factor likely had a significant 
effect, because there was almost no relief when exhumation began and there were only a few 
kilometers of rapid exhumation. As will be shown below, significant relief likely developed after 
the main phase of rapid exhumation had already set the fission-track ages. 

Two of the Shamsi samples, TS167 and TS169, yielded sufficient track-length measurements to 
permit robust thermal modeling (Figure 8); however, given that TS169 yielded twice as many 
measurements, these results are considered more reliable. Both samples were modeled using four 
time-temperature constraints. Models were started at 250 Ma, with a temperature range between 70 
and 160°C, such that all tracks initially formed could be completely annealed. The observed 
shortened track-length distributions require a late-stage heating event, consistent with burial 
heating due to the Oligo-Miocene depositional history discussed above. Three runs were made for 
each sample, with modeled reheating started at 30, 25, or 20 Ma. At this time, the model permitted 
sample temperatures between 25°C and 90°C, consistent with the samples being close to the paleo- 
surface. The peak reheating age of 7 Ma was set to match results from the vertical profile. The 
maximum reheating temperature was constrained to be between 30°C and 100°C. The final 
constraint was the present surface temperature of 1 0-20°C. Acceptable fits were not limited by the 
temperature or the position of the constraints, with the exception of the age of maximum reheating. 
Temperature paths between adjacent constraints had 8 segments. Heating and cooling rates were 
not constrained. Each model run had 10,000 iterations, using a Monte Carlo approach. 

Several conclusions can be drawn from the modeling data (Figure 8 and Table 2). The 
influence of the Miocene sedimentary burial is significant; numerous model runs (not shown) 
which neglect this reheating did not produce good fits. The total annealing temperature (T a ) 
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depends on both apatite kinetic characteristics and cooling rate. T a can be estimated as the best-fit 
temperature at which the oldest preserved track formed; the modeled temperatures can be verified 
by comparison with calculated temperatures reported by Ketcham et al. [1999]. Results from 
shallow samples can be used to understand the behavior of structurally deeper samples with similar 
kinetic characteristics because the temperature where and when the inflection point in the PAZ 
curve formed was approximately T a . This position would have been buried beneath a column of 
bedrock and the Oligo-Miocene sedimentary basin. 

The amount of Miocene reheating experienced by the samples provides a measure of the 
thickness of the Oligocene-Miocene sedimentary basin deposited prior to the onset of uplift of the 
Kyrgyz Range. From this reheating history, the thickness of the exhumed bedrock section and 
hence the subsurface position of the exhumed erosion surface can be estimated. Assuming a 
geothermal gradient of 26°C/km and based on the range of modeled heating of 18-34°C (Table 2) 
between 0.7 and 1.8 km of sediment was deposited at Shamsi. The broad range is partially due to 
the reduced sensitivity of the model at the low temperatures experienced prior to burial. Therefore, 
a more robust calculation taking into account the unsteady thermal conditions within a young 
sedimentary basin is inappropriate. However, the sedimentary thickness calculated is in agreement 
with geological observations from the Chu basin [Chediya, 1 986; Bullen et al., 2003]. 

AFTSolve modeling of sample TS84 from the Boom Gorge cannot precisely constrain the 
timing of final exhumation (Figure 8). The modeling strategy was similar to that used with the 
Shamsi samples, except that reheating began at 1 0, 1 5 or 20 Ma and final cooling began at 2, 3 or 5 
Ma. Best-fit models suggests only 4-22°C of Neogene heating below the sedimentary basin; 
however, acceptable fits permit up to 32°C of heating, consistent with the up to 1.5 km of sediment 
preserved nearby. 

Discussion 
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Multi-compositional age-elevation plots 

Kinetically similar age-elevation curves at Shamsi and Issyk Ata represent an exhumed PA Z 
and are expected to have similar slopes, consistent with geological data that suggests that they have 
experienced similar thermal histories. However, curves within the same transect representing 
apatites with different resistances to annealing are often not parallel because both the T a and the 
influence of Cenozoic burial heating are different. For instance, a more resistant apatite is reset at a 
lower elevation (higher temperature); at intermediate elevations, this apatite will display an older 
age than less resistant apatite because the former has experienced less annealing in the (now 
exhumed) PAZ (Figure 6A). In the case of a single large exhumation event, more resistant apatites 
should lie on a steeper trend within the PAZ. However, subsequent reheating due to burial beneath 
a sedimentary basin can partially reset this exhumed PAZ, creating a zone with a parallel trend 
(Fig. 6C). Subsequent exhumation can expose this history (Fig. 6D). In such a case, the transition 
from parallel to convergent paths in the exhumed PAZ marks the base of older exhumed PAZ. At 
Shamsi, this transition occurs 1 km above the base of the ultimate exhumed PAZ (Fig. 7). Making 
the simplifying assumptions that isotherms have not been perturbed and that the geothermal 
gradient has remained constant, and with independent information on the amount of burial, it is 
possible to estimate the magnitude of the older exhumation event (cf. Fig. 6D). If burial beneath 
the Chu basin was 1 to 1 .5 km, this suggests that the magnitude of the exhumation event associated 
with creation of the Cretaceous erosion surface was 2 to 2.5 km. The range of old ages obtained 
from the Kyrgyz range and the adjacent Chu basin [this study; Bullen et al ., 2001] also implies that 
the penultimate exhumation event was not large enough to reset more thermally resistant all 
apatites. 

Three-stage cooling history 
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Variable age-elevation slopes for kinetically dissimilar apatites support Oligo-Miocene burial 
followed by a three-stage cooling history for the Kyrgyz Range: (1) rapid cooling during removal 
of Cenozoic cover strata, (2) initially reduced erosion and cooling following exhumation of pre- 
Cenozoic bedrock, and (3) renewed rapid cooling and erosion into bedrock (Figure 9). This history 
is best defined by the Shamsi age-elevation curve. Only the lower portion of the Issyk Ata curve 
can presently be well defined with this cooling history. Although erosion has partially removed the 
AFT record of three-stage cooling at Ala Archa, it is also defined by U-Th/He dating of apatite 
there [Bullen et al., 2003]. 

For Issyk Ata and Shamsi, the portion of the age-elevation plot below the inflection point yields 
apparent exhumation rates of ca. 0.4 and ca. 1 km/Myr, respectively. The ca. 1 to 1.5 km of 
sediment overlying the basement would require between 2 to 4 Myr and 1 to 1.5 Myr, respectively, 
to have been removed at these rates. For Shamsi, where the age and paleodepth of the inflection 
point is well constrained as 7 Ma and 4.0 to 4.6 km (Figure 7), the mean exhumation rate from this 
point to the surface is 0.6 to 0.7 km/Myr. This suggests that the rapid cooling recorded in the age- 
elevation profiles reflects the removal of the sedimentary section and that the exhumation rate 
decreased at least slightly afterwards. Similar conclusions about the magnitude of exhumation 
were reached by Bullen et al. [2003] for the Ala Archa section, based on both apatite fission-track 
and U-Th/He thermochronology. In particular, these authors concluded that a ca.l Myr pulse of 
rapid exhumation at ~ 11-10 Ma removed 1.5 km of section at 1.0 to 1.5 km/Myr. Subsequently, 
the exhumation rate decreased to <0.3 km/Myr until 3 Ma, when the rate increased again to 0.8 
km/Myr. These changes were similar in magnitude to synchronous changes in rates of shortening 
along the northern flank of the central Kyrgyz Range and rates of sediment accumulation in the 
Chu Basin [ Bullen et al., 2003] The fission-track data presented in this study cannot address 
whether the exhumation rate changed in the last several Myr at Issyk Ata or Shamsi. 
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The three-stage cooling pattern documented in the Kyrgyz range is likely to be broadly 
applicable. The typical geometry of a foreland-propagating thrust belt includes a foreland basin that 
is eventually disrupted by the advancing reverse-fault-front (Figure 9). In this setting, uplift above a 
new thrust or fold first exhumes the recently deposited sediment of the foreland basin. These 
sediments are typically poorly cemented and therefore easily eroded. Only when the vertical 
displacement of the fault exceeds the thickness of the overlying young sediment are older units 
exposed. Easily eroded, recently deposited sedimentary rocks can be rapidly exhumed with 
negligible surface uplift. As increasing amounts of resistant basement are exposed, exhumation 
rates decrease while surface uplift rates increase, if shortening rates remain steady [c.f., Burbank et 
al., 1999; Sobel and Strecker, 2003]. Surface processes immediately respond to topographic relief 
developed by range growth and initiate drainage networks that can eventually bring erosion rates 
into balance with rock-uplift rates. Orographically enhanced precipitation and glacial erosion may 
further enhance the effectiveness of erosion as the range grows. At this point, if the time scale of 
geomorphic adjustment to range uplift is sufficiently rapid, range morphology may reach a steady 
state dictated by the interplay of surface processes and rock uplift. 

Lateral propagation rate 

The transition from slow cooling, followed by burial, to rapid cooling measures exhumation of 
the Kyrgyz Range via initial erosional stripping of Cenozoic cover strata at the onset of thrusting at 
the range tip. The best estimates for the onset of this rapid exhumation are 1 1 ± 1 Ma at Ala Archa 
[Bullen et al., 2001] and 7 ± 0.5 Ma at Shamsi; these sections are 64 ± 7 km apart (Figure 7). The 
start of exhumation at Issyk Ata, located half-way between these locations, can only be constrained 
as slightly older than 8 Ma. Plotting initiation age versus distance east of Ala Archa shows that the 
range has propagated eastward at ca. 1 6 ± 6 km/Myr from 1 1 to 7 Ma. 
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Progressive changes in the geomorphology of basement rocks that form the easternmost Kyrgyz 
Range provides some independent constraint as to the relative rates of lateral propagation and 
uplift. As soon as rock uplift raises the Cenozoic sedimentary cover rocks above local base level, 
they appear to be rapidly stripped from above the pre-Cenozoic erosion surface and do not 
contribute significantly to the topography of the range. Topographic relief grows primarily within 
the pre-Cenozoic bedrock via surface uplift that results from the competition between rock uplift 
and exhumation. The abrupt increases in mean and maximum basin elevation in the surface uplift 
zone with only nominal gain of internal relief suggests that the form of the nascent range near its 
terminus is dictated primarily by rock uplift above the south-dipping thrust fault, that peak heights 
are correlated with the magnitude of fault slip, and that the magnitude of dissection is minimal. The 
gradient of peak elevations in the surface uplift zone define an east-west slope of 13% ± 2% (Fig. 
4), which, if built by thrusting at a uniform slip rate and without significant erosion of the peaks, 
corresponds to a lateral propagation rate that is only 7 to 9 times the rock uplift rate. Such a 
gradient could not have applied over the entire interval of range growth, however, because it would 
require over 15 km of rock uplift in the central part of the range: an amount inconsistent with the 
Ala Archa fission-track data. Hence, the propagation/rock-uplift ratio appears considerably lower 
near the eastern tip of the range than it does farther west. A quasi-elliptical distribution of fault slip, 
similar to that documented for propagating normal faults [ Dowers , 1993], could be responsible for 
a steeper rock-uplift gradient on the presently propagating easternmost tip of the Kyrgyz Range 
than that derived from Miocene exhumation of the central part of the range. Alternatively, a 
diminished lateral propagation rate could have arose in response to structural interference [Gupta 
and Scholz, 2000 ] with reverse faults and strike-slip faults emerging from the Kungey Alatoo to the 
northeast. Independent geochronological constraints of either fault slip rate or the ages of initiation 
of glaciation in the easternmost Kyrgyz Range as rocks were uplifted above the glacial ELA could 
test these predictions. 
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Geomorphic adjustment to bedrock uplift 

The geomorphic evolution of the Kyrgyz Range supports progressive adjustment of surface 
uplift and incision in response to rock uplift. At the easternmost propagating tip of the Kyrgyz 
Range, limited dissection of the pre-Cenozoic erosion surface enhances surface uplift and dampens 
internal relief. This contrasts with the steady morphology zone in the central part of the range (Fig. 
4) where the pre-Cenozoic erosion surface has been completely removed and large north-facing 
basins with > 1000 m of internal relief predominate. Given the north-vergent underlying thrust 
faults and the structural culmination along the northern edge of the range, the observation that the 
south-facing basins are significantly smaller than the north-facing basins suggests that the northern 
basins have grown at the expense of south-facing ones. This southward migration of the drainage 
divide results in a drainage asymmetry that is opposite to that expected for ranges developed above 
a (singular) actively deforming structure, wherein the steep limb or flank with shorter catchments 
lies above the active fault [ Leeder and Jackson, 1993; Tailing, 1997; Burbank and Anderson, 
2001]. However, this observation is consistent with both range asymmetry driven by > 1000 m 
lower base-level on its north side [Ellis and Densmore, 2003] and with the higher precipitation 
received on the north flank of the range [Aizen et al., 1995]. 

The transition zone, which lies between the zone of surface uplift and the zone of steady 
morphology, provides important insight into the transformation of the Kyrgyz Range by surface 
processes. Peak heights increase gradually from east to west within the transition zone, albeit at a 
much more gentle gradient than in the eastern surface-uplift zone, because of greater competition 
between rock uplift and erosion. The most striking aspect of the transition zone is that mean 
elevation attains a near-constant value of -3200 m despite ongoing growth of peak heights, range 
of hypsometry, and internal relief (Fig. 4). Systematic increase of north-facing basin size explains 
these contradictory elevation and relief trends. The first significant enlargement of north-facing 
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basins occurs where peak elevations exceed 3500 m elevation above mean sea level (Fig. 4), 
enabling Pleistocene glaciation of the range crest to lengthen north-facing basins southward via 
cirque retreat (Figs. 3 and 5), f Oskin and Burbank, submitted manuscript]. Internal relief increases 
only modestly within these glacially expanded north-facing catchments because glacial erosion is 
limited to high elevations. Consequently, hypsometry remains concentrated near the glacial ELA 
[c.f. Brozovic et al., 1997]. Further west within the transition zone, expansion of north-facing 
basins eventually leads to more effective fluvial and glacial erosion at lower elevations, which in 
turn causes the range of hypsometry to expand and internal relief to grow to > 1500m. In the 
westernmost transition zone, the internal relief attains its maximum as rivers in deeply incised 
valleys abut against hillslopes that are at or near the threshold angle for failure. Overall, the 
transition zone appears to be a region of stabilization of elevation, if not locally even of surface 
depression, as large north-facing basins enlarge and mature. 

Coupled Exhumation and Shortening 

Structural reconstructions of the central Kyrgyz Range [Bullen et al., 2003] suggest the pattern 
of shortening may be correlated to the exhumation rate. Initial, rapid shortening corresponds to 
high cooling rates during removal of Cenozoic cover rocks from the crest of the range. This was 
followed by a period of slow shortening rate and backthrusting within the foreland basin, north of 
the Kyrgyz Range, corresponding to a period of slow exhumation of Paleozoic bedrock. In the past 
3 Myr, both shortening rate and exhumation of the bedrock core have increased. If the exhumation 
history of the Kyrgyz Range is driven by the interplay of rock erodability with both rock and 
surface uplift, then this correlation could indicate a feedback mechanism whereby exhumation rates 
are reduced when more resistant rocks are at the surface, thereby building topography. To the 
extent that larger topographic loads require more work for a given increment of shortening and rock 
uplift, shortening and both rock- and surface uplift rates may diminish as topography grows and the 
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locus of deformation shifts elsewhere in the range [e.g., Masek and Duncan, 1998]. If erosion in a 
given range is inefficient and its topography grows, the increased load may cause shortening rates 
to diminish or deformation to shift entirely to a site requiring less work. This model predicts 
variations in exhumation rate with time driven by changes in erosion parameters; the evidence of 
this behavior should be detectable both in the cooling history of the range and in the depositional 
history of the adjacent basin. If the Kyrgyz Range has continued to evolve in a similar manner to 
that expressed in the bedrock and detrital cooling history [Bullen et al. , 2003], then the along-strike 
geomorphic evolution of the range should also be consistent with the proposed exhumation model. 
Here we explore the plausibility of coupling between surface uplift and shortening rate at the range 
scale. 

In the Kyrgyz Range, a short pulse of rapid cooling at several sites is associated with 1 to 1 .5 
km of exhumation; this thickness corresponds well with the thickness of the young sedimentary 
basin which formerly overlay the range. Therefore, most of this pulse of rapid cooling is probably 
attributable to the removal of this sediment during the initial phase of rock uplift. Because the 
young sedimentary rocks are more easily eroded than the underlying Paleozoic units, this short 
episode of rapid cooling likely corresponded to only a small amount of surface uplift [ Bullen et al., 
2003], Only when a significant amount of bedrock was exposed could the surface uplift increase 
and the concurrent exhumation rate decrease. This pattern is seen in the geomorphology of the 
surface uplift zone of the Kyrgyz Range, where easily eroded Cenozoic rocks are rapidly removed 
from the range but the underlying pre-Cenozoic unconformity surface is widely preserved up to 
high elevations. The observation that an initial pulse of rapid exhumation propagates along the 
strike of the range, roughly normal to the principal shortening direction, suggests that the range is 
growing due to lengthening of the bounding structure rather than due to changes in the regional 
stress field. 
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Bullet i et al. [2003] suggested that the decrease in exhumation rate at ca. 1 0 Ma at Ala Archa 
was associated with the observed decrease in range-normal shortening rate and with a shift from 
thrusting to backthrusting. We further suggest that, as the range-bounding thrust fault lengthened, 
this deformation pattern likely propagated along the strike of the range in the same manner. 
Basinward propagating foreland thrust faults and folds, including prominent backthrusts, are 
observed along the length of the range east from Ala Archa. We deduce from the thermal history of 
lateral range propagation that these foreland faults have also propagated eastward as the range itself 
propagated.. Such a spatially varying decrease in shortening rate, if supported from the history of 
these foreland structures, would suggest a feedback between the size of the range and the activity of 
the basal thrust and growing topography [cf. Davis et al. 1983], We speculate that the geomorphic 
transition zone is the surficial expression of a coupled deformation-erosion system that is out of 
balance, leading to a feedback between erosional efficiency and the locus of deformation. 

Steady state time scale 

Consistency between the cooling history and the along-strike geomorphic evolution of the 
Kyrgyz Range supports the hypothesis that exhumation rate, and possibly also shortening rate, are 
modulated by the time scale of adjustment of surface processes to rock uplift. No single steady state 
time scale adequately describes the variety and significance of adjustments toward steady 
morphology in the Kyrgyz Range. Erosion is immediately effective at stripping the cover from the 
pre-Cenozoic erosion surface, supporting an initial, short adjustment time of less than 1 Myr for 
steady state rock uplift and erosion of Cenozoic strata. Conversely, significantly longer time 
periods are necessaiy for erosion to balance rock uplift of resistant Paleozoic basement. In the 
surface uplift zone, below the elevation of glacial ice accumulation, surface processes fail to 
substantively alter the structural form of the range. Low elevation ranges of the western Tien Shan 
commonly preserve extensive areas of the pre-Cenozoic erosion surface [Abdrakmatov et al.. 
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2001], indicating that the steady state time scale in the absence of glaciation is longer than the time 
since Early to Mid Miocene onset of uplift. Glaciation of the easternmost Kyrgyz Range triggers 
the first significant erosional response to rock uplift via establishing, lengthening, and incising 
canyons. Formation of internal relief here corresponds to a sharp decrease in mean basin 
hypsometiy and increasing mean hillslope angles as glacio-fluvial erosion counterbalances, and 
even temporarily exceeds, rock uplift. Over 80% of the adjustment of hypsometry, peak elevations, 
internal relief, and mean slope angles occurs by within the first 25 km of the transition zone, 
suggesting that a quasi-steady state is reached within 2 to 3 Myr after the onset of glaciation and 
canyon cutting in the surface uplift zone. The complete transition to a steady morphology occurs 
gradually and cannot be defined by a sharp break in morphometric indices at the end of the 
transition zone. Mean hypsometry is steady between 3000 and 3500 m for all points west of the 
surface uplift zone. Mean slope angle levels off gradually west of Shamsi, and internal relief 
climbs to maximum values west of Issyk Ata. Overall, approximately 1 10 km of the east Kyrgyz 
Range shows evidence for dynamic adjustment of its geomorphology to balance erosion and rock 
uplift. This probably represents over 10 Myr of range propagation, and suggests that the time scale 
of geomorphic adjustment over active basement-cored uplifts of the Tien Shan spans back into Late 
Miocene time - a significant proportion of the total Cenozoic history of the northern Tien Shan. 

In actively deforming the mountain ranges, the time to steady state is likely to be dependent on 
the efficiency of surface processes and the rates of deformation. Rapid rock-uplift rates tend to 
generate steep hillslopes and rivers, thereby accelerating erosion rates. Numerical models suggest 
that topographic steady state can be attained in <1 Myr where rates of deformation and erosion are 
high (equivalent to several mm/yr: Willett, 1999). Not surprisingly, in the Kyrgyz Range where 
rates are commonly <1 mm/yr, a longer interval is expected to be required to attain steady state. 
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Conclusions 


Comparisons between exhumation histories derived from apatite fission-track analysis and the 
geomorphic evolution based on quantitative DEM analysis of the Kyrgyz Range permits a 
reconstruction of the entire 1 lMyr evolution of the range. Exhumation of the range commenced at 
ca. 1 1 Ma in the vicinity of Ala Archa and has propagated eastward. The lateral propogation rate 
between Ala Archa and Shamsi is ca. 16±6 km/Myr from 1 1 Ma to 7 Ma, and this rate has likely 
slowed as the range propagated farther eastward towards the Kungey Alatoo. Initial rock uplift 
leads to rapid stripping of poorly consolidated, young sediments; rapid exhumation (almost) 
balances rock uplift, creating only limited surface uplift. The ca. 1 .5 km-thick portion of the Chu 
basin which formerly overlay the range above the Cretaceous erosion surface was stripped away in 
1 to 2 Myr, corresponding to the brief episode of rapid cooling documented by apatite fission-track 
analysis and by the zone of surface uplift defined by the presence of the regional erosion surface. 
Once significant amounts of more resistant basement are exposed, the exhumation rate would be 
expected to decrease and surface uplift rate to increase if range-normal shortening rates were 
maintained. However, independent data suggests that the shortening rate actually decreased, 
suggesting that a feedback exists between erosional efficiency and the location and magnitude of 
shortening. In this transition zone, drainage networks propagate into the growing range, eventually 
bringing erosion rates into balance with rock-uplift rates. The timescale for the majority of this 
geomorphic adjustment is ca. 2 to 3 Myr after uplift of the range crest through the glacial ELA and 
development of glacially lengthened and incised canyons, However, additional systematic 
adjustments to mean slope angle, hypsometry, and internal relief occur before the transition to 
steady morphology is complete. Overall, the evolution of the Kyrgyz Range supports a model of a 
systematically time-varying exhumation and shortening rates that are modulated by changes in rock 
erodability, the efficiency of erosion processes, and the time scale of geomorphic adjustment to 
surface uplift. 
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Figures: 


Figure 1. A. Topographic map of central Asia, showing locations of the Tien Shan with respect to 
adjacent plateaus and basins. The Kyrgyz Range forms the northern margin of the central Tien 
Shan. B. Topography of the Kyrgyz Range derived from 3 arcsec digital elevation model [SRTM, 
2003]. Darker shades correspond to lower topography. Locations of apatite fission track (FT) 
samples from Bullen et al. [2001; 2003] and this study shown. Irregular outline defines east half of 
range analyzed in C and in Figure 4. C. Mean slope angle, mean elevation, and maximum 
elevation of the east half of the Kyrgyz Range measured in 1 0 km-wide swaths. Arrows point to 
area of analysis of B contained within white outline. 

Figure 2. Map of the Kyrgyz range and adjacent basins, modified from Mikolaichuk et al. [2003], 
Isopachs show depth to basement; Cenozoic strata of the Chu basin constitutes the majority of this 
sediment. The ranges are composed of Paleozoic units; Mesozoic units are virtually absent. 

Figure 3. Shaded relief image and map of exhumed pre-Cenozoic erosion surface remnants of the 
easternmost Kyrgyz Range. Darker patches are mapped surface remnants with average southward 
dip shown in degrees. Lighter patches show extent of late Pleistocene glaciation shown. Reverse 
fault system bounds northern edge of range. Ages and sample locations from Shamsi River transect 
show complete exhumation of apatite fission-track partial annealing zone northwest of tilted 
erosion surface outcrops. 

Figure 4. A. Topographic characteristics of north-facing basins along the length of the Kyrgyz 
range. Measured basins are shown directly below in B. Peak heights are highest elevations at the 
edge of each basin. Hypsometry shows median elevation bounded by 75 th and 25 th percentile 
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elevations. The difference between these is the hypsometric range, shown shaded as medium grey. 
Internal relief is measured as the difference between the highest and lowest elevations that are the 
same distance upstream from the basin outlet. The 75 th percentile of the distribution of internal 
relief within each basin plotted with area beneath this curve shaded dark grey. Surface uplift zone 
shows sharp increase in peak elevation, hypsometry, and internal relief from east to west in 
proportion to structural growth of the Kyrgyz Range. Adjustment zone shows progressive increases 
in mean slope angle, hypsometric range, and internal relief as north-facing basins expand and incise 
uplifted bedrock. These same morphometric indices are constant in the steady morphology zone. 

Figure 5. Examples of progressive expansion and deepening of north facing basins. Komorchek 
lies at the transition from the surface uplift to the adjustment zone; Tchuk lies within the 
adjustment zone; Ala Archa lies within the steady morphology zone. Erosion from Komorchek to 
Tchuk is dominated by southward expansion of basins, probably via glacial cirque retreat. 
Prominent convexity in stream profile at Tchuk is a result of insufficient fluvial erosion 
downstream of glacially expanded valley. Erosion from Tchuk to Ala Archa dominated by removal 
of convexity by fluvial and glacial incision. Internal relief shown graphically as the shaded region 
between river longitudinal profile and equidistant ridge line elevations, where distance is measured 
up main and tributary streams to the divide. Internal relief effectively captures deepening of basins 
in adjustment zone. 

Figure 6. Schematic illustration of the affect of differing total annealing temperatures (T a ) on 
complex cooling paths. The 4 vertical elevation profiles (A to D) represent a sequence of 
exhumation (cooling) and burial (reheating) events. The upper panel of each pair shows the 
evolving temperature history; the lower panel represents the distribution of apatite fission-track 
ages in a vertical profile at the corresponding time step. Each profile shows the cooling path of 
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apatites that are more (solid line) and less (dashed line) resistant to annealing. The base of the 
partial annealing zone (PAZ) for each type of apatite is considered to be the respective T a . 
Exhumation and burial events are (unrealistically) depicted as instantaneous for clarity. The 
geothermal gradient is assumed to remain constant and advection of isotherms is neglected for 
simplicity. 

A. A large deformation event followed by a long period of quiescence gives shallow samples a 
common age. More resistant apatites are reset at lower elevation (hotter temperature); therefore, 
the slope of the more resistant cooling path is steeper in the PAZ. 

B. Large exhumation event at time tj forms exhumed PAZi. A new PAZ forms at depth. 

C. Section is buried and reheated beneath sediments of a sedimentary basin at time t2. Heating 
resets all apatites below their respective T a . Exhumed PAZ] is partially reset. In the lower part of 
the current PAZ, cooling paths for apatites with different T a may be parallel. If the exhumation 
event was large enough to shift the entire PAZi to a temperature hotter than the appropriate T a , then 
the exhumed PAZI would be completely overprinted. 

D. A second large exhumation event creates new exhumed PAZ2 at time t3. Ongoing exhumation 
could subsequently expose the entire exhumed PAZ in mountainous topography. If the simplifying 
thermal assumptions are correct, the shape of the age-elevation curve can be used to reconstruct the 
magnitude of thermal events. The difference between the magnitude of exhumation ti (a) and 
burial (b) is equal to the difference between the base of PAZi and the base of PAZ2 (c) (for a given 
apatite type). Therefore, if the base of PAZi and PAZ2 can be determined and the magnitude of 
burial can be constrained by thermal modeling or independent geological data, it is possible to 
estimate the magnitude of exhumation that occurred at ti . 
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Figure 7. Fission-track age plotted versus elevation for the Boom Gorge, Shamsi, Issyk Ata, and 
Ala Archa sections (from left to right; east to west). All samples are plotted on a common vertical 
axis. Boom Gorge is in the footwall; other sections are in the hanging wall. Inset in upper right 
shows the young portion of the 3 hanging wall sections. Numbers next to samples indicate kinetic 
character of apatite; italics denote wt%Cl; normal text denotes Dpar. Ala Archa data is from Bullen 
et al. [2003]; Dpar values from this profile are not directly comparable with other Dpar values 
shown as Dpar measurements depend on both the operator and the etching conditions [Sobel et al., 
2004], Upper dashed line indicates the maximum peak elevation. Lower dashed line indicates 
position of base of exhumed PAZ for low-chlorine apatite. Fine dashed line shows approximate 
age-elevation path, constructed following Fig. 6. Amount of exhumation calculated using a 
geothermal gradient of 26°C/km {Gubin, 1986] and the T a (Table 2). Thermal conductivity of 
young sediments and basement were not differentiated. Grey region indicates average thickness of 
sedimentary basin that could have formerly overlain the range; the upper and lower limits are ± ca. 
300 m. The image at the top shows a Multispectral scanner image draped over an SRTM DEM, 
depicting the view of the range looking from the north. White circles denote the location of 
partially reset fission-track samples; grey circles denote fully reset samples. 

Figure 8. Representative track-length models for the Shamsi and Boom Gorge profiles. See text 
for modeling details. Boom Gorge sample modeled using 0.07wt% Cl, chosen by comparison with 
other samples to be equivalent to Dpar of 1.70 pm. 

Figure 9. Schematic evolution of topography in the Kyrgyz range, showing the effect of 
contrasting lithologies on exhumation and surface uplift rates. 

A. Oligocene Early Miocene foreland basin deposited above regional erosion surface. 
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B. Rapid exhumation of sediments in hanging wall of thrust causes rapid cooling but slow surface 
uplift of range. 

C As the area of exposed basement increases relative to young sedimentary cover, erosion rates 
decrease, leading to increasing surface uplift rates. 

D. Range becomes large enough to create an orographic barrier and subsequently place upper 
portion of the range above the Equilibrium Line Altitude (ELA). Glacial erosion causes drainage 
basins on windward side to expand at the expense of more arid, less deeply incised basins on the 
leeward side. 

Figure DR1. Characterization of apatite samples. Compositional data are plotted for all samples. 
Track-length histograms are only plotted for samples with a significant number of measurements. 
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procedures, operating conditions of 15kV beam, beam current of 20 nA (measured on a Faraday cup) and a ca. 15 pm beam diameter. Chlorine was measured for 100 seconds, 
providing a ca. 150 ppm detection limit. Calibration included a Durango apatite standard which was repeatedly measured during the course of analysis. 


Table 2. Summary of AFTSolve model results. 



TS169 

TS167 


min 

max 

min 

max 

Ta (°C) 

112 

113 

100 

103 

Tmin (°C), pre-Cz seds 

35 

50 

57 

57 

Tmax(°C) beneath Cz seds 

81 

84 

75 

82 

burial beneath Cz seds (km) 

1.2 

1.8 

0.7 

1.0 

Wt%Cl 


0.26 


0.08 


Data summarized from all model runs that yielded 
good fits to the observed data. 


